Residue Number System

Question

In the vein of large number challenges I thought this one might be interesting.

In this challenge, we will be using the Residue Number System (RNS) to perform addition, subtraction, and multiplication on large integers.

What is the RNS

The RNS is one of many ways that people have developed to identify integers. In this system, numbers are represented by a sequence of residues (which are the results after a modulus operation (i.e. the remainder after integer division)). In this system, each integer has many representations. To keep things simple, we are going to limit things so that each integer is uniquely represented. I think it is easier to describe what is happening with a concrete example.

Let us look at the first three prime numbers: 2, 3, 5. In the RNS system, we can use these three numbers to uniquely represent any number that is less than 2*3*5 = 30 using residues. Take 21:

21 is less than 30, so we can represent it using the results after modding by 2, 3, and 5. (i.e. the remainder after integer dividing by 2, 3, and 5)

We would identify 21 with the following sequence of integers:

21 ~ { 21 mod 2, 21 mod 3, 21 mod 5 } = {1, 0, 1}

And so in our RNS system, instead of "21", we would use {1,0,1}.

In general given an integer n, we represent n as { n mod 2, ..., n mod p_k } where p_k is the smallest prime such that n is less than the product of all primes less than or equal to p_k.

Another example, say we have 3412. We need to use 2,3,5,7,11,13 here because 2*3*5*7*11*13=30030 whereas, 2*3*5*7*11=2310 which is too small.

3412 ~ { 3412 mod 2, 3412 mod 3, 3412, mod 5, ... , 3412 mod 13 } = { 0, 1, 2, 3, 2, 6 }

You notice that using this system we can represent very large numbers relatively painlessly. Using {1, 2, 3, 4, 5, 6, 7, 8, ... } residues, we can represent numbers up to {2, 6, 30, 210, 2310, 30030, 510510, 9699690 ...} respectively. ( Here is the series )

Our task

We will be using these residues to perform +, -, and * on large numbers. I will describe these processes below. For now here are the input and output specs.

Input

You will be given two (potentially very large) numbers via a stdin or function argument. They will be given as strings of base 10 digits.

For the purposes of outlining the problem further, we call the first input n and the second m. Assume n > m >= 0.

You will also be given + or - or * to indicate the operation to perform.

Output

Let x be an integer. We will use [x] to refer to the RNS representation described above of x.

You are to output [n] <operator> [m] = [result]

How to perform the operations in RNS

These operations are relatively simple. Given two numbers in RNS notation, to add, subtract, or multiply them, simply perform the given operations component-wise and then take the modulus.

i.e.

{ 1, 2, 3 } + { 1, 1, 4 } = { (1+1) mod 2, (2+1) mod 3, (3+4) mod 5 } = {0, 0, 2}

Note that if the number of residues used to represent two different numbers are not the same, when performing operations, you will need to extend the "shorter" number so that it has the same number of residues. This follows the same process. See the test cases for an example.

The same goes if the result requires more residues than either input. Then both inputs need to be "extended".

Important Details

• We will be dealing with big numbers here, but not arbitrarily large. We will be responsible for numbers up to the product of the first 100 primes (see below). To this end, you are given the first 100 primes for free (no byte cost). You may stick them in an array called p or something idiomatic to your language and then subtract the number of bytes used to initiate this array from your final total. This of course means that they may be hard-coded or you may use a built-in to generate them.

• If for someone reason this is the default integer representation used in your language. That is fine.

• You may not use any Arbitrary Precision Integer type unless it is the default of your language. If it is the default, you may not use it to store integers that would not typically fit in 64 bits.

• To be clear, each integer will always be represented with the fewest residues possible. This goes for both input and output.

• I think the other specs should prevent this, but to be redundant: you may not perform the given operation on the inputs and then and then change everything to RNS and then output. You must change the inputs to RNS and then perform the operations to produce the output.

Test Cases

1. Input:

n = 10
m = 4
+

Output:

{ 0, 1, 0 } + { 0, 1 } = { 0, 2, 4 }

Explanation:

First, change each number to its RNS representation as described above:

10 ~ {0,1,0} and 4 ~ {0,1}. Notice that when we want to do component-wise addition, that 10 has more components than 4. Therefore we must "extend" the shorter number. So we will briefly write 4 ~ {0,1} --> {0,1, 4 mod 5} = {0,1,4}. Now we proceed with addition and then take the modulus.

2. Input

n=28
m=18
+

Output:

 [ 0, 1, 3 ] + [0, 0, 3 ] = [ 0, 1, 1, 4 ]

3. Input (me mashing my face on the keyboard)

n=1231725471982371298419823012819231982571923
m=1288488183
*

Output (broken onto separate lines for readability):

[1, 2, 3, 6, 2, 10, 2, 1, 12, 16, 7, 15, 34, 29, 31, 5, 55, 32, 66, 61, 3, 76, 52, 14, 65, 44, 99, 57 ] 
* 
[1, 0, 3, 3, 4, 8, 9, 10, 8, 0 ] 
= 
[1, 0, 4, 4, 8, 2, 1, 10, 4, 0, 17, 7, 27, 21, 44, 51, 56, 9, 6, 9, 12, 0, 52, 36, 43, 68, 99, 24, 96, 39, 96, 66, 125] 

n requires 28 primes. m requires 10. n*m requires 33.

4. Input

n=8709668761379269784034173446876636639594408083936553641753483991897255703964943107588335040121154680170867105541177741204814011615930342030904704147856733048115934632145172739949220591246493529224396454328521288726490
m=1699412683745170450115957274739962577420086093042490863793456500767137147999161679589295549397604032154933975242548831536518655879433595016
-

Output:

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 509]
-
[0, 2, 1, 6, 1, 12, 11, 18, 14, 28, 21, 36, 37, 42, 16, 52, 41, 60, 16, 70, 49, 78, 80, 88, 49, 100, 13, 106, 4, 112, 68, 130, 36, 138, 37, 150, 0, 162, 8, 172, 163, 180, 18, 192, 129, 198, 135, 222, 78, 228, 90, 238, 57, 250, 36, 262, 87, 270, 206, 280, 193, 292, 253, 310, 224, 316, 57, 336, 48, 348]
=
[0, 1, 4, 1, 10, 1, 6, 1, 9, 1, 10, 1, 4, 1, 31, 1, 18, 1, 51, 1, 24, 1, 3, 1, 48, 1, 90, 1, 105, 1, 59, 1, 101, 1, 112, 1, 0, 1, 159, 1, 16, 1, 173, 1, 68, 1, 76, 1, 149, 1, 143, 1, 184, 1, 221, 1, 182, 1, 71, 1, 90, 1, 54, 1, 89, 1, 274, 1, 299, 1, 266, 1, 228, 1, 340, 1, 170, 1, 107, 1, 340, 1, 88, 1, 157, 1, 143, 1, 22, 1, 22, 1, 58, 1, 296, 1, 371, 1, 140]

n uses 100 primes. m uses 70 primes. n-m uses 99 primes.

I checked these using the ChineseRem built-in implementation of the Chinese Remainder theorem on GAP (which basically takes RNS numbers and changes them to base 10 integers). I believe they are correct. If something seems fishy, please let me know.

---

For those who care, the product of the first 100 primes is:

471193079990618495316248783476026042202057477340967552018863483961641533584503
422120528925670554468197243910409777715799180438028421831503871944494399049257
9030720635990538452312528339864352999310398481791730017201031090

This number is 1 larger than the maximal number we can represent using the given system (and 100 prime limitation).

Somewhat related

Answer 1

GAP

Some Background: I will admit that when I created this question, so many months ago, I did not have a method for solving the hard part of this question: determining the correct number of primes to use. We have a lot of very intelligent people on this site, and I really expected that someone would figure out a way to do it fairly quickly. However, as this did not happen, I wasn't even sure that it was truly possible to solve this problem. So, I had to take the time to devise a method. I believe that what I have done does not break any of the rules in this challenge, of course I would love this to be fact checked.

I also slightly regret the choice of code-golf because the solutions are slightly more in depth than would usually fit the tag format. Having said this, to follow site rules, there is a "golfed" version of my solution at the bottom of this post.

---

Code

### The first 100 primes;
primes := Primes{[1..100]};

### In many of the functions below, the 'string' variable is a string of digits
###


### Returns the 'index' digit of 'string' as an integer
GetValueAsInt := function(string, index) 
    return IntChar(string[index]) - 48;
end;

### Used in the 'modulus' function. See that comment for more information. 
### Calculates the contribution to the modulus of a digit 'digit' in the 10^power place.
### 'integer' is the modulus
digit_contribution := function(digit, integer, power)
    local result, i;
    result := 1;
    for i in [0..power-1] do
        result := ( result * (10 mod integer) ) mod integer;
    od;
    result := (result * (digit mod integer) ) mod integer;
    return result;
end;

### This modulus function is used to calculate the modulus of large numbers without storing them
##### as large numbers.
### It does so by breaking them into digits, and calculating the contribution of each digit.
### e.g. 1234 mod 5 = (1000 mod 5)(1 mod 5) + (200 mod 5)(2 mod 5) + (10 mod 5)(3 mod 5) + (4 mod 5)
### It actually mods after every calculation to ensure that we never get a number larger
##### than the modulus ('integer') squared, which will never be even close to 10^64-1
modulus := function(string, integer)
    local i, result, digit, len;
    len := Length(string);
    result := 0;
    for i in [1..len] do
        digit :=  IntChar(string[i]) -48;
        result := ( result + digit_contribution(digit, integer, len-i) )  mod integer;
    od;
    return result;
end;

### This returns the product of the first i-1 primes (mod j). It must not (and does not)
##### ever store an integer larger than 2^64-1
phi_i := function(i,j)
    local index, result;
    result := 1;
    for index in [1..i-1] do
        result := ( result * primes[index] ) mod primes[j];
    od;
    return result;
end;

### Calculates the first residues of 'string' mod the first 100 primes
get_residues := function(string) 
    local p, result;
    result := [];
    for p in primes do
        Add( result, modulus(string, p) );  
    od; 
    return result;
end;

### Gets the ith element in the partial_chinese array, given the previous elements
### See the explanation section and partial_chinese function for more info
get_partial_i := function( i, residues, previous_array )
    local index, result;
    result := residues[i];
    for index in [1..Length(previous_array)] do
        result := ( result - previous_array[index]*phi_i(index,i) ) mod primes[i]; 
    od;     
    result := ( result / phi_i(i,i) ) mod primes[i];
    return result;
end;

### returns an array such that the sum of prod_primes(i)*array[i] is equal to the integer value
##### that is represented by the residues. (It basically just does the CRT without
##### actually summing everything.) prod_primes(i) is the product of the first i-1 primes 
### See the explanation for a bit more info
### This is what allows us to determine the minimal number of primes to represent a RNS number
partial_chinese := function( string )
    local array, i, residues;
    residues := get_residues(string);
    array := [];        
    for i in [1 .. Length(primes)] do
        Add( array, get_partial_i( i, residues, array ) );
    od;
    return array;   
end;

### Same as partial_chinese but takes input in a different form.
partial_chinese_from_residues := function(residues)
    local array, i;
    array := [];        
    for i in [1 .. Length(primes)] do
        Add( array, get_partial_i( i, residues, array ) );
    od;
    return array;
end;

### gives you the number of primes needed to represent an integer. Basically asks how 
##### many trailing zeros there are in the chinese array.
get_size := function(string)
    local array, i, len, result;
    array := partial_chinese(string);
    len := Length(array);
    for i in [0..len-1] do
        if  not (array[len-i] = 0) then
            return len -i;
        fi; 
    od; 
    Print("ERROR: get_size().\n");
    return 0;
end;

### Same as above but with different input format
get_size_from_residues := function(residues)
    local array, i, len, result;
    array := partial_chinese_from_residues(residues);
    len := Length(array);
    for i in [0..len-1] do
        if  not (array[len-i] = 0) then
            return len -i;
        fi; 
    od; 
    Print("ERROR: get_size().\n");
    return 0;
end;

### the actual function. inputs are all strings
f := function(in1, in2, opperation)
    local residues_1, residues_2, residues_result, i;
    residues_1 := get_residues(in1);
    residues_2 := get_residues(in2);
    residues_result := [];
    if opperation = "+" then
        for i in [1..Length(primes)] do
            Add( residues_result, ( residues_1[i] + residues_2[i] ) mod primes[i]);
        od;     
    elif opperation = "*" then
        for i in [1..Length(primes)] do
            Add( residues_result, ( residues_1[i] * residues_2[i] ) mod primes[i]);
        od;     
    elif opperation = "-" then
        for i in [1..Length(primes)] do
            Add( residues_result, ( residues_1[i] - residues_2[i] ) mod primes[i]);
        od;     
    fi;
    Print(residues_1{[1..get_size(in1)]}, " ", opperation, " ", residues_2{[1..get_size(in2)]}, " = ", residues_result{[1..get_size_from_residues(residues_result)]} );
end;

Explanation:

To start off, we calculate all the 100 of the residues for both inputs. We do this with the modulus function in the code. I tried to be careful so that we use the built-in mod function after every step. This ensures that we never have a number that is larger than 540^2, which is 1 less than the 100th prime squared.

After we have all the residues, we can perform the given operation and mod each entry again. Now we have a unique designator for the result, but we need to determine the minimal number of entries that we have to use to represent the result and each of the inputs.

Figuring out how many residues we actually need is by far the hardest part of this problem. To determine this, we perform most of the steps of the Chinese Remainder Theorem (CRT). However, we obviously have to make modifications so that we don't end up with numbers that are too large.

Let prod(i) be the sum of the first i-1 primes. For example,

prod(1) = 1
prod(2) = 2
prod(3) = 6
prod(4) = 30
etc

Let X be an integer. Let {r_i} be the residues of X, that is

r_i = X mod p_i

Where p_i is the ith prime. This is for 1<i<=100 in our case.

Now we are going to use the CRT to find a sequence {u_i} such that the sum over i of prod(i) * u_i is equal to X. Note that each u_i is also technically a residue, as u_i < p_i. Furthermore, if X < prod(i) then u_i = 0. This is of critical importance. It means that by examining trailing zeros, we can determine how many of the residues r_i we actually need to represent X in the RNS.

If you care to examine some sequences of u_i, the partial_chinese function returns the u_i sequence.

By messing around with the CRT, I was able to find a recursive formula for the u_i values, solving the issue of determining how many residues we need.

The formula is:

u_i = [ r_i - SUM ] / prod(i)       (mod p_i)

Where SUM is the sum over j in [1,i) of u_j * prod(i).

Of course, prod(i) cannot actually be calculated in some cases because it is too large. For this purpose, I used the phi_i function. This function returns prod(j) (mod p_i). It mods at every step, so we never actually calculate anything that is too large.

If you are curious where this formula comes from, I would recommend working a couple of CRT examples, which can be found on the wikipedia page.

Finally, for each input as well as our output, we calculate the u_i sequence and then determine the trailing zeros. Then we throw out that many r_i from the end of the residue sequences.

---

"Golfed" Code, 2621 bytes

primes:=Primes{[1..100]};GetValueAsInt:=function(string,index)return IntChar(string[index])-48;end;digit_contribution := function(digit, integer, power)local result, i;result:=1;for i in [0..power-1] do result := ( result * (10 mod integer) ) mod integer;od;result:=(result*(digit mod integer) ) mod integer;return result;end;modulus:=function(string, integer)local i,result,digit,len;len:=Length(string);result:=0;for i in [1..len] do digit:= IntChar(string[i])-48;result:=(result+digit_contribution(digit,integer,len-i)) mod integer;od;return result;end;phi_i:=function(i,j)local index,result;result:=1;for index in [1..i-1] do result:=(result*primes[index] ) mod primes[j];od;return result;end;get_residues:=function(string) local p,result;result:=[];for p in primes do Add(result,modulus(string,p));od;return result;end;get_partial_i:=function(i,residues,previous_array)local index,result;result:=residues[i];for index in [1..Length(previous_array)] do result:=(result-previous_array[index]*phi_i(index,i) ) mod primes[i];od;result:=(result/phi_i(i,i)) mod primes[i];return result;end;partial_chinese:=function(string)local array,i,residues;residues:=get_residues(string);array:=[];for i in [1 .. Length(primes)] do Add(array,get_partial_i(i,residues,array));od;return array;end;partial_chinese_from_residues:=function(residues)local array,i;array:=[];for i in [1..Length(primes)] do Add(array,get_partial_i(i,residues,array));od;return array;end;get_size:=function(string)local array,i,len,result;array:=partial_chinese(string);len:=Length(array);for i in [0..len-1] do if not (array[len-i] = 0) then return len-i;fi;od;Print("ERROR: get_size().\n");return 0;end;get_size_from_residues:=function(residues)local array,i,len,result;array:=partial_chinese_from_residues(residues);len:=Length(array);for i in [0..len-1] do if not (array[len-i]=0) then return len-i;fi;od;Print("ERROR: get_size().\n");return 0;end;f:=function(in1,in2,opperation)local residues_1,residues_2,residues_result,i;residues_1:=get_residues(in1);residues_2:=get_residues(in2);residues_result:=[];if opperation = "+" then for i in [1..Length(primes)] do Add(residues_result,(residues_1[i]+residues_2[i] ) mod primes[i]);od;elif opperation = "*" then for i in [1..Length(primes)] do Add(residues_result,(residues_1[i]*residues_2[i])mod primes[i]);od;elif opperation = "-" then for i in [1..Length(primes)] do Add(residues_result,(residues_1[i]-residues_2[i]) mod primes[i]);od;fi;Print(residues_1{[1..get_size(in1)]}, " ", opperation, " ", residues_2{[1..get_size(in2)]}, " = ", residues_result{[1..get_size_from_residues(residues_result)]} );end;

Answer 2

Mathematica, not golfed

rns[d_,l_]:=Table[Reap[
    FoldPairList[Sow@QuotientRemainder[10#+#2,Prime@i]&,0,d]
  ][[2,1,-1,2]],{i,l}];
plus[a_,b_]:=Mod[a+b,Prime@Range@Length@a];
subtract[a_,b_]:=Mod[a-b,Prime@Range@Length@a];
times[a_,b_]:=Mod[a b,Prime@Range@Length@a];
mag[f_]:=LengthWhile[FoldList[#/#2&,f,Prime@Range@100],#>1.1&];
ext[m_,n_,i_]:=Fold[Mod[1##,Prime@i]&,m,Prime@Range@n];
multi[e_,p_,t_]:=Tr@Position[Mod[e Range@p,p],p-t];
appx[d_] := N@FromDigits[{d~Take~UpTo[6], Length@d}]

• Function rns[d_,l_] converts a base-10 integer d into a RNS integer of length l.

• Function plus / times / subtract add / multiply / subtract one RNS integer to / from another, both of which are of the same length.

• Function mag[f_] estimates the approximate magnitude of the floating point number f in terms of the lower bound of the length of its RNS representation.

• Function ext[m_,n_,i_] finds out the remainder from division of the product of m and Prime[Range@n] by Prime[i].

• Function multi[e_,p_,t_] gives the smallest multiplier m satisfying that Divisible[m*e+t,p]

• Function appx[d_] take the first 6 digits of a decimal integer and gives its approximate floating point value.

---

With the help of functions above, now we are able to solve a tricky problem - to determine the length of the result.

First, I have to clarify that it is not an easy task to determine the RNS length of an integer. For small integers, we can compare them directly with the product of primes. But for very large integers, since it is forbidden to calculate the product of primes infinitely accurate, such comparison no longer works.

For example, given that the product of prime 1 to 30 is 3.16*10^46, the RNS length of integers around 3.16*10^46 can possibly be 29 or 30. The function mag will give 29 as a reference for these integers, showing that both 29 and 30 are possible.

Once knowing the magnitude, we simply represent the integer according to that magnitude directly, hoping to calculate its true length. The trick here is to add some new numbers to the original number and modify its RNS representation, until the representation is all zero.

For example, mag[211.] is 4, and its length 4 representation is {1, 1, 1, 1}.

step 1:   {1,1,1,1} -> {0,2,2,2}  by adding  (1) * 1 = 1
step 2:   {0,2,2,2} -> {0,0,1,6}  by adding  (2) * 2 = 4
step 3:   {0,0,1,6} -> {0,0,0,2}  by adding  (2*3) * 4 = 24
step 4:   {0,0,0,2} -> {0,0,0,0}  by adding  (2*3*5) * 6 = 180
step 5:   calculate 211 + (1 + 4 + 24 + 180) ~ 420

By adding some number, we increase 211 to the smallest number that is divisible by 210 (2*3*5*7). And now we conclude that the original number is greater than 210, since 420 equals "approximately" two times of 210. It is not difficult to imagine that if we start from 209, the final number is "approximately" 210.

Function length[f_,n_] takes the floating point value f to estimate the magnitude, and correct it base on its RNS representation n.

length[f_,n_]:=With[{g=mag@f},
    g+If[#==0,1,Round[(#+f)/Times@@Prime@Range@g]-1]&[
      FoldList[Times,1.,Prime[Range[g-1]]].
      FoldPairList[
        Block[{i=#2,m},
          {m=multi[ext[1,i-1,i],Prime@i,Part@##],rnsPlus[#,ext[m,i-1,#]&/@Range[g]]}
        ]&,n,Range[g]]]]

---

Function rnsOperation[a_,b_,op_,rnsop_] gives rnsop[a,b] and op is corresponding normal operations used to get approximate results based on floating point values.

rnsOperation[a_,b_,op_,rnsop_]:=Block[{c=op[appx@a,appx@b],m},
    m=mag[c];m=length[c,rnsop[rns[a,m],rns[b,m]]];rnsop[rns[a,m],rns[b,m]]]

Example

rnsOperation[
    IntegerDigits@1231725471982371298419823012819231982571923,
    IntegerDigits@1288488183,
    Times, times]
(* {1,0,4,4,8,2,1,10,4,0,17,7,27,21,44,51,56,9,6,9,12,0,52,36,43,68,99,24,96,39,96,66,125} *)

Answer 3

Python 3, 435 bytes

This challenge has been on my bucket list for a while, but it's only recently that: a) I put time and attention into actually attempting an answer; and b) actually tested my idea to calculate the size of the numbers (and so the number of primes by comparing it to the size of the primorials) using some unholy combination of logarithms and the Chinese Remainder Theorem. Unfortunately, working with logarithms while trying to determine the number of primes that, for example, large_primorial + 3 requires, meant that I was having to find ways to work around floating-point problems.

And so, this is a port of Liam's answer.

Try it online!

from functools import reduce as R
G=range
d=lambda s:[R(lambda z,c:(z*10+int(c))%q,s,0)for q in p]
h=lambda j,i:R(lambda z,q:z*q%p[i],p[:j],1)
def s(r):
 a=[];z=99
 for i in G(100):
  P=p[i];u=r[i]
  for j in G(len(a)):u=(u-a[j]*h(j,i))%P
  for k in G(1,P):
   if h(i,i)*k%P<2:break
  a+=u*k%P,
 while(a[z]<1)*z:z-=1
 return r[:z+1]
def f(a,b,n):u=d(a);v=d(b);print(s(u),n,s(v),'=',s([eval(str(u[i])+n+str(v[i]))%p[i]for i in G(100)]))

Explanation

While trying to port Liam's answer, I personally found some of the explanation given was confusingly worded, so this is my attempt at explaining his algorithm.

First, we get the residues of n and m.

res1 = get_residues(n)
res2 = get_residues(m)

This involves turning all the digits of the input strings and turning them into numbers modulo each of our primes, e.g. for 28, we would have [(20 + 8) mod 2, (20 + 8) mod 3, (20 + 8) mod 5, etc]

def get_residues(string):
    result = []
    for p in primes:
        result.append(reduce(lambda z, c:(z*10+int(c)) % p, string, 0))

Then we add, multiply, or subtract the residues pairwise using eval()

result = []
for i in range(len(primes)):
    result.append((eval(str(res1[i]) + op + str(res2[i])) % primes[i])

Then we get the sizes of our residues, i.e. the minimum number of primes we need.

size1 = get_size(res1)
size2 = get_size(res2)
size3 = get_size(result)

Getting the sizes is the trickiest and most code-intensive part. We use the partial_chinese function, which gets us our sequence of u_i to determine the size with. More on u_i in a moment.

def get_size(residues):
    array = partial_chinese(residues)
    size = len(residues)-1
    while array[size] == 0 and size:
        size -= 1
    return size+1  # to prevent off-by-one errors from 0-indexing

The sequence u_i is calculated by taking each residue r_i, subtracting the sum u_j * primorial(j) for j in [1, i), and then dividing by primorial(i), all modulo primes[i]. That is, u_i = (r_i - SUM) / primorial(i). More on our primorial and division functions in a moment.

def partial_chinese(residues):
    array = []
    for i in range(len(primes)):
        array.append(get_partial_i(i, residues, array))
    return array

def get_partial_i(i, residues, previous_array):
    result = residues[i]
    for j in range(len(previous_array)):
        result = (result - previous_array[j] * phi_i(j, i)) % primes[i]
    result = result * inverse(phi_i(i, i), primes[i]) % primes[i]
    return result

phi_i(j, i) calculates primorial(j) mod primes[i]. Division modulo any prime p is easily implemented by checking for multiplicative inverses manually, as we can be sure that any possible u_i is 0 <= u_i < p is guaranteed to be coprime to p and so is guaranteed a multiplicative inverse.

def phi_i(j, i):
    return reduce(lambda z, q: z * q % primes[i], primes[:j], 1)

def inverse(n, p):
    for i in range(1, p):
        if n * i % p == 1:
            return i

With all that done, we print our string and we're done.

print(res1[:size1], op, res2[:size2], "=", result[:size3])

What's next

This was fun to implement. I still want to see if I can use logarithms in some way in another answer. And I'd like to implement this code or something like in a functional golfing language, like APL or Jelly. Any and all golfing suggestions and corrections are welcome!