A limited set of basic number theoretic tools

Various tools for working with integers and their factors, primes, congruences, etc.

Over the years, I've put together a few useful tools for working with integers, congruences, factors of numbers, divisors, etc. None of these tools requires the use of the symbolic toolbox, although I've tried where possible to make them accessible to any numeric class, including doubles, syms, int64 and my own VPIJ tools.

In this suite of tools, you will find functions to list all factor pairs of a number, the modular inverse of a number, a rrefgf tool, Legendre and Jacobi symbols, discrete logarithms, modular square roots, the Euler totient function, etc.

For example, modSqrt solves the problem:

mod(x^2,n) = a

taken modulo 685, there are two solutions to this problem.

x = modSqrt(15,685)
x =
   120   565
   
mod(x.^2,685)
ans =
    15    15

modInv solves the problem

mod(a*x,n) = 1

x = modInv(12,137)
x =
    80

mod(x*12,137)
ans =
     1

There are 6 pairs of numbers that are factor pairs of the number 36 (if we assume the first element of the pair is never larger than the second.)

factorpairs(36)
ans =
     1    36
     2    18
     3    12
     4     9
     6     6

modLog solves what is often called the discrete logarithm problem. It finds a solution to the problem

mod(a^x,n) == b

Note that this solution need not be unique.

modLog(11,760,9931)
ans =
   137

So, we would have

mod(11^137,9931) == 760

Of course you cannot compute that result using mod, but powerMod will suffice.

powerMod(11,137,9931)
ans =
   760

rrefgf perfroms row manipulations on a matrix in the ring of integers under multiplication modulo some number N (sometimes called GF(N).)

A = randi(2,[4,4]) - 1
A =
     1     1     0     1
     0     0     1     1
     0     0     0     1
     1     0     1     1
[Arref,Abasis,Arank,Ainv] = rrefgf(A,2)
Arref =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
Abasis =
     1     2     3     4
Arank =
     4
Ainv =
     0     1     0     1
     1     1     1     1
     0     1     1     0
     0     0     1     0

So this random matrix has full rank in GF(2). As a test of the inverse produced:

mod(A*Ainv,2)
ans =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

There are tools for working with quadratic residues. In this example, P happens to be prime.

P = 10000019;

legendreSymbol(12345,P)
ans =
    -1
    
legendreSymbol(12346,P)
ans =
     1

If the legendreSymbol (a/p) is -1, then a is not a quadratic residue.

a = 12345;
legendreSymbol(a,P)
ans =
    -1

In that case, no modular square root will exist, as a is not a quadratic residue.

modSqrt(a,P)
Warning: No solution for the congruential sqrt exists 
In modSqrt (line 123) 
ans =
     []

However, 12346 IS a quadratic residue.

legendreSymbol(a+1,P)
ans =
     1
modSqrt(a+1,P)
ans =
     4722543     5277476

You can find a list of all quadratic residues of a number. (This list will be long for bigger numbers, so be careful.)

quadraticResidues(100)
ans =
  Columns 1 through 15
     0     1     4     9    16    21    24    25    29    36    41    44    49    56    61
  Columns 16 through 22
    64    69    76    81    84    89    96

The Jacobi Symbol (a/P) works just like the Legendre symbol, but for composite moduli.

Another simple tool is properDivisors, which is here only because it is considerably faster than the similar symbolic toolbox tool divisors.

properDivisors(210)
ans =
     1     2     3     5     6     7    10    14    15    21    30    35    42    70   105

Another cute tool is the function repetend. It computes the repeating part of any rational number, composed as the ratios ot two integers. It can even handle different bases.

For example, the ratio 1/107 is a repeating decimal in base 10. It looks like this, but the repeating pattern is much longer then a double precision number can identify.

1/107
ans =
       0.00934579439252336

[Rstr,Nrstr] = repetend(1,107)
Rstr =
    '00934579439252336448598130841121495327102803738317757'
Nrstr =
  1×0 empty char array

There are 53 decimal digits in the repeat pattern.

numel(Rstr)
ans =
    53

To verify that claim, we can use syms.

vpa(sym(1)/107,150)
ans =
0.00934579439252336448598130841121495327102803738317757009345794392523364485981308411214953271028037383177570093457943925233644859813084112149532710280374

In base 12, the fractions 1/12, 1/6, 1/3, and 1/2 all are finitely exactly represented with a non-repeating (so entirely zero extension) radix. So 1/6, in base would look like 0.200000000...

[Rstr,Nrstr] = repetend(1,6,12)
Rstr =
  0×0 empty char array
Nrstr =
    '2'

Enjoy.