[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](https://en.wikipedia.org/wiki/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](https://oeis.org/A002110) )
**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](https://codegolf.stackexchange.com/questions/71730/arbitrary-precision-integer-division)
  
[Answer]
# 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](/questions/tagged/code-golf "show questions tagged '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 `i`th 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 `mod`s 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](https://en.wikipedia.org/wiki/Chinese_remainder_theorem).
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.
---
[Answer]
# 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]
# [Python 3](https://docs.python.org/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](https://codegolf.stackexchange.com/a/91392/47581).
[Try it online!](https://tio.run/##XVTLciJJDLzzFXVx0A3tiHqoXu3p81wdcyU4tA0s@NHgBjwx/Lw3s@wNMxthqiSVlEolFoc/p@1@cB@HbmEb5RrlGxUbZQw@8Azt3CgL2@J2iDvEBLcgJrA9q/AWEAvwo2kmSqmIcEQ44U64M6E0cTWBdXEzW9G1dF3pWiw@CA//iWY8w4Gp7GFiIceExKJcDsZKm0zKnMEW4gS3tozhCpolvuUQ1peDCQS3gQ@RscgYwW0qw@NwJO2I60jaUR3CuSILiTspB5WiLI66OBJ2JOxI2BHOURFHrkJJhEIIhCaa2CJuOSgxuUoRm7hCrkKthYSF4EJcSbQohHB8rz9H9YT2hPTUwotZfmzG/avanIfH037/clS718N@PKlxvTo/rlV/VL8mP7uxH/5ZT1bdS//6sOrVsV38qr7sS/PYVpeZ0fPdcKoe6/rmrTk2ut7sR/WmdoM6LCfb/wqfml17VfnWXmZvN4fFbtkcFu3TsjH1ZLXeqGM11u1E9d1ieXfpcp4oou2I9rMyWvNR3XcsvDt3Iy74THn6THlZD1Vf1@25q863/eJpOdtWaA1u91@Jz19YzX3BUruN2lY7pMyeb@5/2PZhXPfPeOjn3ZkhqPd7u3tZV/3isvxh6tmlvdx2ZgKdTudxUOOivczNsrDfVH3z0AxsvwKNu3dcD/XdYaRAx@pcN0NzrN7rZtpNYSzW7/1LdTyN1RmD1PNhTvuddl20@d/sCH9sqqmZNuUzn9aTK/f2b3f27eorl7ZcFVvY8bvaJpanq4TMwBWgLQ2uIS0WwHqJJifrorE5STG1wd5kvMLxkRYrbUqSkknuGyFFnUNIMWDpM1YvJtFOsNsijAaHv@yziE46ueyC9y7g3TuBm03K6O@jdjlIxsbo6FNyzmssFdbfA0WbqFOIRnv855sYoxirJRlkmGB8dmhotdNZS9RiJCYfonNM4St2zBrxmDOioWRrtc/GCvo5b@EKe3txNmHHMGK0eKJIJiDf2JBcZL0Wr4nogSTECpAmojXYabAQi7IUoIMTH7zWEaQhKlYbgwb86GWfss0YI@N3I2BEMKMDz1vB94DRDfTxJkGnFCEIfn7QNXx@yx//Ag "Python 3 – 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!
]