////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 1994-2021 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or .
//
// This file is part of Octave.
//
// Octave is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// Octave is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Octave; see the file COPYING. If not, see
// .
//
////////////////////////////////////////////////////////////////////////
#if defined (HAVE_CONFIG_H)
# include "config.h"
#endif
#include
#include
#include
#include
#include
#include
#include
#include "lo-ieee.h"
#include "mx-base.h"
#include "oct-base64.h"
#include "oct-binmap.h"
#include "oct-time.h"
#include "quit.h"
#include "Cell.h"
#include "data.h"
#include "defun.h"
#include "error.h"
#include "errwarn.h"
#include "interpreter-private.h"
#include "oct-map.h"
#include "ov-class.h"
#include "ov-complex.h"
#include "ov-cx-mat.h"
#include "ov-cx-sparse.h"
#include "ov-float.h"
#include "ov-flt-complex.h"
#include "ov-flt-cx-mat.h"
#include "ov.h"
#include "ovl.h"
#include "pager.h"
#include "parse.h"
#include "pt-mat.h"
#include "utils.h"
#include "variables.h"
#include "xnorm.h"
OCTAVE_NAMESPACE_BEGIN
// This function implements the Schrage technique for modular multiplication.
// The returned value is equivalent to "(a*b) % modulus"
// but calculated without overflow.
uint64_t safemultiply (uint64_t a, uint64_t b, uint64_t modulus)
{
uint64_t q = modulus / a;
uint64_t r = modulus - q*a;
uint64_t term1 = a * (b % q);
uint64_t term2 = (r < q) ? r * (b/q) : safemultiply (r, b / q, modulus);
return (term1 > term2) ? (term1 - term2) : (term1 + modulus - term2);
}
// This function returns "pow(a,b) % modulus"
// but calculated without overflow.
uint64_t safepower (uint64_t a, uint64_t b, uint64_t modulus)
{
uint64_t retval = 1;
while (b > 0)
{
if (b & 1)
retval = safemultiply (retval, a, modulus);
b >>= 1;
a = safemultiply (a, a, modulus);
}
return retval;
}
// This function implements a single round of Miller-Rabin primality testing.
// Returns false if composite, true if pseudoprime for this divisor.
bool millerrabin (uint64_t div, uint64_t d, uint64_t r, uint64_t n)
{
uint64_t x = safepower (div, d, n);
if (x == 1 || x == n-1)
return true;
for (uint64_t j = 1; j < r; j++)
{
x = safemultiply (x, x, n);
if (x == n-1)
return true;
}
return false;
}
DEFUN (__isprimelarge__, args, ,
doc: /* -*- texinfo -*-
@deftypefn {} {@var{x} =} __isprimelarge__ (@var{n})
Use the Miller-Rabin test to find out whether the scalar N is prime or
composite. The input N is required to be a scalar 64-bit integer. You
probably want to call isprime(N) instead of directly calling this function.
@seealso{isprime, factor}
@end deftypefn */)
{
int nargin = args.length ();
if (nargin != 1)
print_usage ();
// This function is intended for internal use by isprime.m,
// so the following error handling should not be necessary. But it is
// probably good practice for any curious users calling it directly.
uint64_t n = args(0).xuint64_scalar_value
("__isprimelarge__: unable to convert input into uint64_t scalar. Call isprime() instead.");
if (n <= 37) // n is too small to test reliably with this code.
// This case happens only if the user is directly calling this function.
// Calling it from isprime.m will prevent this case before it gets here.
error ("__isprimelarge__: input argument too small. Call isprime() instead.");
// rewrite n as d * 2^r + 1, for odd d
uint64_t d = n-1;
uint64_t r = 0;
while ( ! (d & 1))
{
d >>= 1;
r++;
}
// Miller-Rabin test with the first 12 primes.
// If the number passes all 12 tests, then it is prime.
// If it fails any, then it is composite.
// The first 12 primes suffice to test all 64-bit integers.
if ( ! millerrabin ( 2, d, r, n)) return ovl (false);
if ( ! millerrabin ( 3, d, r, n)) return ovl (false);
if ( ! millerrabin ( 5, d, r, n)) return ovl (false);
if ( ! millerrabin ( 7, d, r, n)) return ovl (false);
if ( ! millerrabin (11, d, r, n)) return ovl (false);
if ( ! millerrabin (13, d, r, n)) return ovl (false);
if ( ! millerrabin (17, d, r, n)) return ovl (false);
if ( ! millerrabin (19, d, r, n)) return ovl (false);
if ( ! millerrabin (23, d, r, n)) return ovl (false);
if ( ! millerrabin (29, d, r, n)) return ovl (false);
if ( ! millerrabin (31, d, r, n)) return ovl (false);
if ( ! millerrabin (37, d, r, n)) return ovl (false);
return ovl (true); // if all the way here, then it is prime
/*
Mathematical references for the curious as to why we need only
the 12 smallest primes for testing all 64-bit numbers:
(1) https:oeis.org/A014233
Comment: a(12) > 2^64. Hence the primality of numbers < 2^64 can be
determined by asserting strong pseudoprimality to all prime bases <= 37
(=prime(12)). Testing to prime bases <=31 does not suffice,
as a(11) < 2^64 and a(11) is a strong pseudoprime
to all prime bases <= 31 (=prime(11)). - Joerg Arndt, Jul 04 2012
(2) https:arxiv.org/abs/1509.00864
Strong Pseudoprimes to Twelve Prime Bases
Jonathan P. Sorenson, Jonathan Webster
In addition, a source listed here: https://miller-rabin.appspot.com/
reports that all 64-bit numbers can be covered with only 7 divisors,
namely 2, 325, 9375, 28178, 450775, 9780504, and 1795265022.
I couldn't find a peer-reviewed article to back it up though,
so I am using the 12 primes <= 37 in this code.
*/
}
/*
%!assert (__isprimelarge__ (uint64 (12345)), false)
%!assert (__isprimelarge__ (uint64 (2147483647)), true)
%!assert (__isprimelarge__ (uint64 (2305843009213693951)), true)
%!assert (__isprimelarge__ (uint64 (18446744073709551557)), true)
*/
OCTAVE_NAMESPACE_END