Solovay-Strassen method of Primality Test

Try it on GfG Practice

We have already been introduced to primality testing in the previous articles in this series. 

• Introduction to Primality Test and School Method
• Fermat Method of Primality Test
• Primality Test | Set 3 (Miller–Rabin)

The Solovay–Strassen test is a probabilistic algorithm used to check if a number is prime or composite (not prime). It is not 100% accurate but gives a high probability of correctness when run multiple times. Before jumping into the test, let's break down some key concepts:

To perform the test, we need to understand two important mathematical symbols:

• Legendre Symbol (a/p)
• Jacobian Symbol (a/n)

These symbols help us determine properties of numbers in a structured way.

Legendre Symbol (a/p)

This is a special way to compare two numbers:

• a: Any integer
• p: A prime number

It is written as (a/p) and tells us one of the following three things:

1. (a/p) = 0 → If a is divisible by p, meaning a % p == 0
2. (a/p) = 1 → if there exists an integer k such that k² = a(mod p)
3. (a/p) = -1 → Otherwise (if the above two cases are false)

Example

• If p = 7, and a = 9, then:
  (9/7) = 1 because 9 has a square root 3 (since 3² ≡ 9).
• If a = 2, then (2/7) = -1 because 2 is not a square modulo 7.

Mathematician Euler discovered a useful shortcut to compute (a/p):

(a/p) = a^(p-1)/2 mod p

This formula helps in quickly checking the Legendre symbol.

Jacobian Symbol (a/n)

The Jacobian symbol is a generalization of the Legendre symbol, but it works when n is not necessarily a prime. It is calculated by breaking down n into prime factors and applying the Legendre symbol to each factor.

If n is a prime number, then:

(a/n) = (a/p)

So, Jacobian symbol = Legendre symbol when n is prime.

The Solovay–Strassen Primality Test

Now that we understand these symbols, let's go step by step through the primality test.

Step 1: Choose a Random Number

• Pick a random number a, such that 2 ≤ a ≤ n - 1.
• This number a acts as a "witness" to help determine whether n is prime.

Step 2: Check the GCD (Greatest Common Divisor)

• Compute gcd(a, n) (the largest number that divides both a and n).
• If gcd(a, n) > 1, then n is definitely composite (not prime).
• This step ensures n has no common factors with a.

Step 3: Compute Two Values

We compute two important values:

1. a^((n-1)/2) mod n (using fast exponentiation)
2. (a/n) (using the Jacobian symbol)

Step 4: Compare the Two Values

• If these two values do not match, then n is composite.
• If they are equal, then n is probably prime.

Why Do We Repeat This Test Multiple Times?

Since this is a probabilistic test, we cannot be 100% sure that n is prime after one test. So, we repeat the test with different random values of a multiple times.

• If the test fails even once, n is definitely composite.
• If the test passes many times, then n is probably prime.

// C++ program to implement Solovay-Strassen Primality Test
#include <bits/stdc++.h>
using namespace std;

// Function to perform modular exponentiation
long long modulo(long long base, long long exponent, long long mod) {
    long long x = 1, y = base;

    while (exponent > 0) {
        if (exponent % 2 == 1) {
            x = (x * y) % mod;
        }
        y = (y * y) % mod;
        exponent /= 2;
    }
    return x % mod;
}

// Function to calculate the Jacobian symbol (a/n)
int calculateJacobian(long long a, long long n) {
    if (a == 0) {
        return 0; // (0/n) = 0
    }

    int ans = 1;

    if (a < 0) {
        a = -a; // (a/n) = (-a/n) * (-1/n)
        if (n % 4 == 3) {
            ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
        }
    }

    if (a == 1) {
        return ans; // (1/n) = 1
    }

    while (a) {
        if (a < 0) {
            a = -a; // (a/n) = (-a/n) * (-1/n)
            if (n % 4 == 3) {
                ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
            }
        }

        while (a % 2 == 0) {
            a /= 2;
            if (n % 8 == 3 || n % 8 == 5) {
                ans = -ans;
            }
        }

        swap(a, n);

        if (a % 4 == 3 && n % 4 == 3) {
            ans = -ans;
        }

        a %= n;

        if (a > n / 2) {
            a -= n;
        }
    }

    return (n == 1) ? ans : 0;
}

// Function to perform the Solovay-Strassen Primality Test
bool solovoyStrassen(long long p, int iterations) {
    if (p < 2 || (p % 2 == 0 && p != 2)) {
        return false;
    }

    for (int i = 0; i < iterations; i++) {
        long long a = rand() % (p - 1) + 1;
        long long jacobian = (p + calculateJacobian(a, p)) % p;
        long long mod = modulo(a, (p - 1) / 2, p);

        if (!jacobian || mod != jacobian) {
            return false;
        }
    }
    return true;
}

int main() {
    int iterations = 50;
    long long num1 = 15, num2 = 13;

    if (solovoyStrassen(num1, iterations)) {
        printf("%lld is prime\n", num1);
    } else {
        printf("%lld is composite\n", num1);
    }

    if (solovoyStrassen(num2, iterations)) {
        printf("%lld is prime\n", num2);
    } else {
        printf("%lld is composite\n", num2);
    }

    return 0;
}

// Java program to implement Solovay-Strassen Primality Test

import java.util.Random;

class GFG {
    
    // Modulo function to perform binary exponentiation
    static long modulo(long base, long exponent, long mod) {
        long x = 1;
        long y = base;

        while (exponent > 0) {
            if (exponent % 2 == 1) {
                x = (x * y) % mod;
            }
            y = (y * y) % mod;
            exponent /= 2;
        }
        return x % mod;
    }

    // Function to calculate the Jacobian symbol of a given number
    static long calculateJacobian(long a, long n) {
        if (n <= 0 || n % 2 == 0) {
            return 0;
        }

        long ans = 1;
        if (a < 0) {
            a = -a; // (a/n) = (-a/n) * (-1/n)
            if (n % 4 == 3) {
                ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
            }
        }
        if (a == 1) {
            return ans; // (1/n) = 1
        }

        while (a != 0) {
            if (a < 0) {
                a = -a; // (a/n) = (-a/n) * (-1/n)
                if (n % 4 == 3) {
                    ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
                }
            }
            while (a % 2 == 0) {
                a /= 2;
                if (n % 8 == 3 || n % 8 == 5) {
                    ans = -ans;
                }
            }

            long temp = a;
            a = n;
            n = temp;

            if (a % 4 == 3 && n % 4 == 3) {
                ans = -ans;
            }
            a %= n;
            if (a > n / 2) {
                a -= n;
            }
        }
        return (n == 1) ? ans : 0;
    }

    // Function to perform the Solovay-Strassen Primality Test
    static boolean solovayStrassen(long p, int iteration) {
        if (p < 2) {
            return false;
        }
        if (p != 2 && p % 2 == 0) {
            return false;
        }

        Random rand = new Random();
        for (int i = 0; i < iteration; i++) {
            long r = Math.abs(rand.nextLong());
            long a = r % (p - 1) + 1;
            long jacobian = (p + calculateJacobian(a, p)) % p;
            long mod = modulo(a, (p - 1) / 2, p);

            if (jacobian == 0 || mod != jacobian) {
                return false;
            }
        }
        return true;
    }

    // Driver code
    public static void main(String[] args) {
        int iter = 50;
        long num1 = 15;
        long num2 = 13;

        if (solovayStrassen(num1, iter)) {
            System.out.println(num1 + " is prime");
        } else {
            System.out.println(num1 + " is composite");
        }

        if (solovayStrassen(num2, iter)) {
            System.out.println(num2 + " is prime");
        } else {
            System.out.println(num2 + " is composite");
        }
    }
}

# Python3 program to implement Solovay-Strassen 
# Primality Test 
import random

# modulo function to perform binary exponentiation
def modulo(base, exponent, mod):
    x = 1
    y = base
    while exponent > 0:
        if exponent % 2 == 1:
            x = (x * y) % mod

        y = (y * y) % mod
        exponent //= 2

    return x % mod

# To calculate Jacobian symbol of a given number
def calculateJacobian(a, n):
    if a == 0:
        return 0  # (0/n) = 0

    ans = 1
    if a < 0:
        # (a/n) = (-a/n)*(-1/n)
        a = -a
        if n % 4 == 3:
            # (-1/n) = -1 if n = 3 (mod 4)
            ans = -ans

    if a == 1:
        return ans  # (1/n) = 1

    while a:
        if a < 0:
            # (a/n) = (-a/n)*(-1/n)
            a = -a
            if n % 4 == 3:
                # (-1/n) = -1 if n = 3 (mod 4)
                ans = -ans

        while a % 2 == 0:
            a //= 2
            if n % 8 == 3 or n % 8 == 5:
                ans = -ans

        # Swap
        a, n = n, a

        if a % 4 == 3 and n % 4 == 3:
            ans = -ans
        a %= n

        if a > n // 2:
            a -= n

    if n == 1:
        return ans

    return 0

# To perform the Solovay-Strassen Primality Test
def solovoyStrassen(p, iterations):
    if p < 2:
        return False
    if p != 2 and p % 2 == 0:
        return False

    for _ in range(iterations):
        # Generate a random number a
        a = random.randrange(1, p)
        jacobian = (p + calculateJacobian(a, p)) % p
        mod = modulo(a, (p - 1) // 2, p)

        if jacobian == 0 or mod != jacobian:
            return False

    return True

if __name__ == "__main__":
    iterations = 50
    num1 = 15
    num2 = 13

    if solovoyStrassen(num1, iterations):
        print(num1, "is prime")
    else:
        print(num1, "is composite")

    if solovoyStrassen(num2, iterations):
        print(num2, "is prime")
    else:
        print(num2, "is composite")

// C# program to implement the Solovay-Strassen Primality Test

using System;

class GfG {
    // Function to perform modular exponentiation
    static long Modulo(long Base, long exponent, long mod) {
        long x = 1;
        long y = Base;

        while (exponent > 0) {
            if (exponent % 2 == 1) {
                x = (x * y) % mod;
            }
            y = (y * y) % mod;
            exponent /= 2;
        }
        return x % mod;
    }

    // Function to calculate the Jacobian symbol of a given number
    static long CalculateJacobian(long a, long n) {
        if (n <= 0 || n % 2 == 0) {
            return 0;
        }

        long ans = 1;
        if (a < 0) {
            a = -a; // (a/n) = (-a/n) * (-1/n)
            if (n % 4 == 3) {
                ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
            }
        }
        if (a == 1) {
            return ans; // (1/n) = 1
        }

        while (a != 0) {
            if (a < 0) {
                a = -a; // (a/n) = (-a/n) * (-1/n)
                if (n % 4 == 3) {
                    ans = -ans; // (-1/n) = -1 if n ≡ 3 (mod 4)
                }
            }
            while (a % 2 == 0) {
                a /= 2;
                if (n % 8 == 3 || n % 8 == 5) {
                    ans = -ans;
                }
            }

            long temp = a;
            a = n;
            n = temp;

            if (a % 4 == 3 && n % 4 == 3) {
                ans = -ans;
            }
            a %= n;
            if (a > n / 2) {
                a -= n;
            }
        }
        return (n == 1) ? ans : 0;
    }

    // Function to perform the Solovay-Strassen Primality Test
    static bool SolovayStrassen(long p, int iteration) {
        if (p < 2) {
            return false;
        }
        if (p != 2 && p % 2 == 0) {
            return false;
        }

        Random rand = new Random();
        for (int i = 0; i < iteration; i++) {
            long r = Math.Abs(rand.Next());
            long a = r % (p - 1) + 1;
            long jacobian = (p + CalculateJacobian(a, p)) % p;
            long mod = Modulo(a, (p - 1) / 2, p);

            if (jacobian == 0 || mod != jacobian) {
                return false;
            }
        }
        return true;
    }

    // Driver Code
    static void Main() {
        int iter = 50;
        long num1 = 15;
        long num2 = 13;

        if (SolovayStrassen(num1, iter)) {
            Console.WriteLine(num1 + " is prime");
        } else {
            Console.WriteLine(num1 + " is composite");
        }

        if (SolovayStrassen(num2, iter)) {
            Console.WriteLine(num2 + " is prime");
        } else {
            Console.WriteLine(num2 + " is composite");
        }
    }
}

// JavaScript program to implement Solovay-Strassen Primality Test

// Modulo function to perform binary exponentiation
function modulo(base, exponent, mod) {
    let x = 1;
    let y = base;
    
    while (exponent > 0) {
        if (exponent % 2 == 1) {
            x = (x * y) % mod;
        }
        y = (y * y) % mod;
        exponent = Math.floor(exponent / 2);
    }
    return x % mod;
}

// To calculate Jacobian symbol of a given number
function calculateJacobian(a, n) {
    if (n <= 0 || n % 2 == 0) {
        return 0;
    }
    
    let ans = 1;
    if (a < 0) {
        a = -a;
        if (n % 4 == 3) {
            ans = -ans;
        }
    }
    
    if (a == 1) {
        return ans;
    }
    
    while (a !== 0) {
        if (a < 0) {
            a = -a;
            if (n % 4 == 3) {
                ans = -ans;
            }
        }
        
        while (a % 2 == 0) {
            a = Math.floor(a / 2);
            if (n % 8 == 3 || n % 8 == 5) {
                ans = -ans;
            }
        }
        
        let temp = a;
        a = n;
        n = temp;
        
        if (a % 4 == 3 && n % 4 == 3) {
            ans = -ans;
        }
        
        a %= n;
        if (a > Math.floor(n / 2)) {
            a -= n;
        }
    }
    
    return n === 1 ? ans : 0;
}

// To perform the Solovay-Strassen Primality Test
function solovoyStrassen(p, iteration) {
    if (p < 2 || (p % 2 === 0 && p !== 2)) {
        return false;
    }
    
    for (let i = 0; i < iteration; i++) {
        let a = Math.floor(Math.random() * (p - 2)) + 1;
        let jacobian = (p + calculateJacobian(a, p)) % p;
        let mod = modulo(a, Math.floor((p - 1) / 2), p);
        
        if (jacobian === 0 || mod !== jacobian) {
            return false;
        }
    }
    return true;
}

// Driver Code
let iter = 50;
let num1 = 15;
let num2 = 13;

console.log(num1 + (solovoyStrassen(num1, iter) ? " is prime" : " is composite"));
console.log(num2 + (solovoyStrassen(num2, iter) ? " is prime" : " is composite"));

Output : 

15 is composite
13 is prime

Time Complexity: Using fast algorithms for modular exponentiation, the running time of this algorithm is O(k*(log n)³), where k is the number of different values we test.
Auxiliary Space: O(1) as it is using constant space for variables

Accuracy: It is possible for the algorithm to return an incorrect answer. If the input n is indeed prime, then the output will always probably be correctly prime. However, if the input n is composite, then it is possible for the output to probably be incorrect prime. The number n is then called an Euler-Jacobi pseudoprime.