Code for 大数取模 modulo arithmetic

2024年5月19日


/// <summary>
/// Find b^p % m.
/// </summary>
/// <param name="b"></param>
/// <param name="p"></param>
/// <param name="m"></param>
/// <returns></returns>
public static int FindModularExponentiation(int b, int p, int m)
{
	int p1 = (int)Math.Ceiling(Math.Log(m) / Math.Log(b));
	if (p1 <= p)
		logger.LogInformation($"{b}^{p1} >= {m}. ({p1} <= {p})");
	else
		return (int)NumberService.PowerWithSwitch(b, p);

	int q = Math.DivRem(p, p1, out int r);
	logger.LogInformation($"{b}^{p} mod {m} = {b}^{{{p1}*{q}+{r}}} mod {m} = {b}^{{{p1}*{q}}} * {b}^{r} mod {m} ");

	Debug.Assert(NumberService.PowerWithSwitch(b, r) < m);
	logger.LogInformation($"= ({b}^{p1} mod {m})^{q} * {b}^{r} mod {m} ");

	// Even thought p1 is miminal power that makes b^p1 >= m,
	// b^p1 can still be multiples of m, if b is large.
	int res = (int)(NumberService.PowerWithSwitch(b, p1) % m);
	logger.LogInformation($"= {res}^{q} * {b}^{r} mod {m} ");


	Dictionary<int, int> factors = NumberService.Factor(res);
	SortedDictionary<int, int> sortedFactors = new SortedDictionary<int, int>(factors.ToDictionary(kv => kv.Key, kv => kv.Value * q));

	if (r > 0)
	{
		if (sortedFactors.TryGetValue(b, out int c))
			sortedFactors[b] = c + r;
		else
			sortedFactors[b] = r;
	}

	logger.LogInformation("= {0} mod {1}",
		string.Join(" * ", sortedFactors.Select(kv => kv.Value == 1 ? kv.Key.ToString() : kv.Key + "^" + kv.Value)),
		m);

	logger.LogInformation("Then we divide the expression in {0} parts.", sortedFactors.Count);


	List<int> remainders = new List<int>(sortedFactors.Count);
	int result = 1;
	foreach (var f in sortedFactors)
	{
		int remainder;
		using (logger.BeginScope(" "))
		{
			remainder = FindModularExponentiation(f.Key, f.Value, m);
		}
		remainders.Add(remainder);
		result = (int)((long)result * remainder % m);
	}
	logger.LogInformation(string.Join(" * ", remainders) + " mod " + m + " = " + result);
	logger.LogInformation("\n");
	return result;
}

There are a couple of places to watch out.

To start with, you can’t use Math.Pow() because it works on double and double is inherently inaccurate. ^[1] You can use BigInteger.Pow() or use a switch-based solution^[2] (the return value must be long).

Next, as I talked about in the previous article, you have to do prime factorization. There is certainly room to optimize the factorization.

The following is the log.

 2^30 >= 1000000007. (30 <= 3500)
 2^3500 mod 1000000007 = 2^{30*116+20} mod 1000000007 = 2^{30*116} * 2^20 mod 1000000007
 = (2^30 mod 1000000007)^116 * 2^20 mod 1000000007
 = 73741817^116 * 2^20 mod 1000000007
 = 2^20 * 101^116 * 491^116 * 1487^116 mod 1000000007
 Then we divide the expression in 4 parts.
  101^5 >= 1000000007. (5 <= 116)
  101^116 mod 1000000007 = 101^{5*23+1} mod 1000000007 = 101^{5*23} * 101^1 mod 1000000007
  = (101^5 mod 1000000007)^23 * 101^1 mod 1000000007
  = 510100431^23 * 101^1 mod 1000000007
  = 3^23 * 101 * 617^23 * 275581^23 mod 1000000007
  Then we divide the expression in 4 parts.
 |  3^19 >= 1000000007. (19 <= 23)
 |  3^23 mod 1000000007 = 3^{19*1+4} mod 1000000007 = 3^{19*1} * 3^4 mod 1000000007
 |  = (3^19 mod 1000000007)^1 * 3^4 mod 1000000007
 |  = 162261460^1 * 3^4 mod 1000000007
 |  = 2^2 * 3^4 * 5 * 251 * 32323 mod 1000000007
 |  Then we divide the expression in 5 parts.
 |  4 * 81 * 5 * 251 * 32323 mod 1000000007 = 143178169
 |

 |  617^4 >= 1000000007. (4 <= 23)
 |  617^23 mod 1000000007 = 617^{4*5+3} mod 1000000007 = 617^{4*5} * 617^3 mod 1000000007
 |  = (617^4 mod 1000000007)^5 * 617^3 mod 1000000007
 |  = 924113713^5 * 617^3 mod 1000000007
 |  = 617^3 * 4127^5 * 223919^5 mod 1000000007
 |  Then we divide the expression in 3 parts.
 | |  4127^3 >= 1000000007. (3 <= 5)
 | |  4127^5 mod 1000000007 = 4127^{3*1+2} mod 1000000007 = 4127^{3*1} * 4127^2 mod 1000000007
 | |  = (4127^3 mod 1000000007)^1 * 4127^2 mod 1000000007
 | |  = 291595893^1 * 4127^2 mod 1000000007
 | |  = 3 * 4127^2 * 97198631 mod 1000000007
 | |  Then we divide the expression in 3 parts.
 | |  3 * 17032129 * 97198631 mod 1000000007 = 830680711
 | |

 | |  223919^2 >= 1000000007. (2 <= 5)
 | |  223919^5 mod 1000000007 = 223919^{2*2+1} mod 1000000007 = 223919^{2*2} * 223919^1 mod 1000000007
 | |  = (223919^2 mod 1000000007)^2 * 223919^1 mod 1000000007
 | |  = 139718211^2 * 223919^1 mod 1000000007
 | |  = 3^2 * 1693^2 * 27509^2 * 223919 mod 1000000007
 | |  Then we divide the expression in 4 parts.
 | |  9 * 2866249 * 756745081 * 223919 mod 1000000007 = 649279648
 | |

 |  234885113 * 830680711 * 649279648 mod 1000000007 = 773529357
 |

 |  275581^2 >= 1000000007. (2 <= 23)
 |  275581^23 mod 1000000007 = 275581^{2*11+1} mod 1000000007 = 275581^{2*11} * 275581^1 mod 1000000007
 |  = (275581^2 mod 1000000007)^11 * 275581^1 mod 1000000007
 |  = 944887036^11 * 275581^1 mod 1000000007
 |  = 2^22 * 275581 * 236221759^11 mod 1000000007
 |  Then we divide the expression in 3 parts.
 | |  236221759^2 >= 1000000007. (2 <= 11)
 | |  236221759^11 mod 1000000007 = 236221759^{2*5+1} mod 1000000007 = 236221759^{2*5} * 236221759^1 mod 1000000007
 | |  = (236221759^2 mod 1000000007)^5 * 236221759^1 mod 1000000007
 | |  = 34449048^5 * 236221759^1 mod 1000000007
 | |  = 2^15 * 3^10 * 478459^5 * 236221759 mod 1000000007
 | |  Then we divide the expression in 4 parts.
 | | |  478459^2 >= 1000000007. (2 <= 5)
 | | |  478459^5 mod 1000000007 = 478459^{2*2+1} mod 1000000007 = 478459^{2*2} * 478459^1 mod 1000000007
 | | |  = (478459^2 mod 1000000007)^2 * 478459^1 mod 1000000007
 | | |  = 923013085^2 * 478459^1 mod 1000000007
 | | |  = 5^2 * 907^2 * 203531^2 * 478459 mod 1000000007
 | | |  Then we divide the expression in 4 parts.
 | | | |  203531^2 >= 1000000007. (2 <= 2)
 | | | |  203531^2 mod 1000000007 = 203531^{2*1+0} mod 1000000007 = 203531^{2*1} * 203531^0 mod 1000000007
 | | | |  = (203531^2 mod 1000000007)^1 * 203531^0 mod 1000000007
 | | | |  = 424867674^1 * 203531^0 mod 1000000007
 | | | |  = 2 * 3 * 7 * 11 * 311 * 2957 mod 1000000007
 | | | |  Then we divide the expression in 6 parts.
 | | | |  2 * 3 * 7 * 11 * 311 * 2957 mod 1000000007 = 424867674
 | | | |

 | | |  25 * 822649 * 424867674 * 478459 mod 1000000007 = 549840858
 | | |

 | |  32768 * 59049 * 549840858 * 236221759 mod 1000000007 = 476999544
 | |

 |  4194304 * 275581 * 476999544 mod 1000000007 = 256423819
 |

  143178169 * 101 * 773529357 * 256423819 mod 1000000007 = 390957170


  491^4 >= 1000000007. (4 <= 116)
  491^116 mod 1000000007 = 491^{4*29+0} mod 1000000007 = 491^{4*29} * 491^0 mod 1000000007
  = (491^4 mod 1000000007)^29 * 491^0 mod 1000000007
  = 120048155^29 * 491^0 mod 1000000007
  = 5^29 * 23^29 * 1043897^29 mod 1000000007
  Then we divide the expression in 3 parts.
 |  5^13 >= 1000000007. (13 <= 29)
 |  5^29 mod 1000000007 = 5^{13*2+3} mod 1000000007 = 5^{13*2} * 5^3 mod 1000000007
 |  = (5^13 mod 1000000007)^2 * 5^3 mod 1000000007
 |  = 220703118^2 * 5^3 mod 1000000007
 |  = 2^2 * 3^2 * 5^3 * 36783853^2 mod 1000000007
 |  Then we divide the expression in 4 parts.
 | |  36783853^2 >= 1000000007. (2 <= 2)
 | |  36783853^2 mod 1000000007 = 36783853^{2*1+0} mod 1000000007 = 36783853^{2*1} * 36783853^0 mod 1000000007
 | |  = (36783853^2 mod 1000000007)^1 * 36783853^0 mod 1000000007
 | |  = 832054252^1 * 36783853^0 mod 1000000007
 | |  = 2^2 * 208013563 mod 1000000007
 | |  Then we divide the expression in 2 parts.
 | |  4 * 208013563 mod 1000000007 = 832054252
 | |

 |  4 * 9 * 125 * 832054252 mod 1000000007 = 244107792
 |

 |  23^7 >= 1000000007. (7 <= 29)
 |  23^29 mod 1000000007 = 23^{7*4+1} mod 1000000007 = 23^{7*4} * 23^1 mod 1000000007
 |  = (23^7 mod 1000000007)^4 * 23^1 mod 1000000007
 |  = 404825426^4 * 23^1 mod 1000000007
 |  = 2^4 * 23 * 97^4 * 1409^4 * 1481^4 mod 1000000007
 |  Then we divide the expression in 5 parts.
 | |  1409^3 >= 1000000007. (3 <= 4)
 | |  1409^4 mod 1000000007 = 1409^{3*1+1} mod 1000000007 = 1409^{3*1} * 1409^1 mod 1000000007
 | |  = (1409^3 mod 1000000007)^1 * 1409^1 mod 1000000007
 | |  = 797260915^1 * 1409^1 mod 1000000007
 | |  = 5 * 11 * 127 * 157 * 727 * 1409 mod 1000000007
 | |  Then we divide the expression in 6 parts.
 | |  5 * 11 * 127 * 157 * 727 * 1409 mod 1000000007 = 340621374
 | |

 | |  1481^3 >= 1000000007. (3 <= 4)
 | |  1481^4 mod 1000000007 = 1481^{3*1+1} mod 1000000007 = 1481^{3*1} * 1481^1 mod 1000000007
 | |  = (1481^3 mod 1000000007)^1 * 1481^1 mod 1000000007
 | |  = 248367620^1 * 1481^1 mod 1000000007
 | |  = 2^2 * 5 * 17 * 19 * 1481 * 38447 mod 1000000007
 | |  Then we divide the expression in 6 parts.
 | |  4 * 5 * 17 * 19 * 1481 * 38447 mod 1000000007 = 832442651
 | |

 |  16 * 23 * 88529281 * 340621374 * 832442651 mod 1000000007 = 453927012
 |

 |  1043897^2 >= 1000000007. (2 <= 29)
 |  1043897^29 mod 1000000007 = 1043897^{2*14+1} mod 1000000007 = 1043897^{2*14} * 1043897^1 mod 1000000007
 |  = (1043897^2 mod 1000000007)^14 * 1043897^1 mod 1000000007
 |  = 720938986^14 * 1043897^1 mod 1000000007
 |  = 2^14 * 1043897 * 360469493^14 mod 1000000007
 |  Then we divide the expression in 3 parts.
 | |  360469493^2 >= 1000000007. (2 <= 14)
 | |  360469493^14 mod 1000000007 = 360469493^{2*7+0} mod 1000000007 = 360469493^{2*7} * 360469493^0 mod 1000000007
 | |  = (360469493^2 mod 1000000007)^7 * 360469493^0 mod 1000000007
 | |  = 474109271^7 * 360469493^0 mod 1000000007
 | |  = 43^7 * 1297^7 * 8501^7 mod 1000000007
 | |  Then we divide the expression in 3 parts.
 | | |  43^6 >= 1000000007. (6 <= 7)
 | | |  43^7 mod 1000000007 = 43^{6*1+1} mod 1000000007 = 43^{6*1} * 43^1 mod 1000000007
 | | |  = (43^6 mod 1000000007)^1 * 43^1 mod 1000000007
 | | |  = 321363007^1 * 43^1 mod 1000000007
 | | |  = 7 * 29 * 43 * 1153 * 1373 mod 1000000007
 | | |  Then we divide the expression in 5 parts.
 | | |  7 * 29 * 43 * 1153 * 1373 mod 1000000007 = 818609210
 | | |

 | | |  1297^3 >= 1000000007. (3 <= 7)
 | | |  1297^7 mod 1000000007 = 1297^{3*2+1} mod 1000000007 = 1297^{3*2} * 1297^1 mod 1000000007
 | | |  = (1297^3 mod 1000000007)^2 * 1297^1 mod 1000000007
 | | |  = 181825059^2 * 1297^1 mod 1000000007
 | | |  = 3^2 * 13^2 * 1297 * 4662181^2 mod 1000000007
 | | |  Then we divide the expression in 4 parts.
 | | | |  4662181^2 >= 1000000007. (2 <= 2)
 | | | |  4662181^2 mod 1000000007 = 4662181^{2*1+0} mod 1000000007 = 4662181^{2*1} * 4662181^0 mod 1000000007
 | | | |  = (4662181^2 mod 1000000007)^1 * 4662181^0 mod 1000000007
 | | | |  = 931524616^1 * 4662181^0 mod 1000000007
 | | | |  = 2^3 * 11 * 101 * 311 * 337 mod 1000000007
 | | | |  Then we divide the expression in 5 parts.
 | | | |  8 * 11 * 101 * 311 * 337 mod 1000000007 = 931524616
 | | | |

 | | |  9 * 169 * 1297 * 931524616 mod 1000000007 = 63530421
 | | |

 | | |  8501^3 >= 1000000007. (3 <= 7)
 | | |  8501^7 mod 1000000007 = 8501^{3*2+1} mod 1000000007 = 8501^{3*2} * 8501^1 mod 1000000007
 | | |  = (8501^3 mod 1000000007)^2 * 8501^1 mod 1000000007
 | | |  = 341771203^2 * 8501^1 mod 1000000007
 | | |  = 41^2 * 8501 * 8335883^2 mod 1000000007
 | | |  Then we divide the expression in 3 parts.
 | | | |  8335883^2 >= 1000000007. (2 <= 2)
 | | | |  8335883^2 mod 1000000007 = 8335883^{2*1+0} mod 1000000007 = 8335883^{2*1} * 8335883^0 mod 1000000007
 | | | |  = (8335883^2 mod 1000000007)^1 * 8335883^0 mod 1000000007
 | | | |  = 944903287^1 * 8335883^0 mod 1000000007
 | | | |  = 71 * 97 * 137201 mod 1000000007
 | | | |  Then we divide the expression in 3 parts.
 | | | |  71 * 97 * 137201 mod 1000000007 = 944903287
 | | | |

 | | |  1681 * 8501 * 944903287 mod 1000000007 = 904205081
 | | |

 | |  818609210 * 63530421 * 904205081 mod 1000000007 = 524796728
 | |

 |  16384 * 1043897 * 524796728 mod 1000000007 = 982403768
 |

  244107792 * 453927012 * 982403768 mod 1000000007 = 580316551


  1487^3 >= 1000000007. (3 <= 116)
  1487^116 mod 1000000007 = 1487^{3*38+2} mod 1000000007 = 1487^{3*38} * 1487^2 mod 1000000007
  = (1487^3 mod 1000000007)^38 * 1487^2 mod 1000000007
  = 288008282^38 * 1487^2 mod 1000000007
  = 2^38 * 1487^2 * 144004141^38 mod 1000000007
  Then we divide the expression in 3 parts.
 |  2^30 >= 1000000007. (30 <= 38)
 |  2^38 mod 1000000007 = 2^{30*1+8} mod 1000000007 = 2^{30*1} * 2^8 mod 1000000007
 |  = (2^30 mod 1000000007)^1 * 2^8 mod 1000000007
 |  = 73741817^1 * 2^8 mod 1000000007
 |  = 2^8 * 101 * 491 * 1487 mod 1000000007
 |  Then we divide the expression in 4 parts.
 |  256 * 101 * 491 * 1487 mod 1000000007 = 877905026
 |

 |  144004141^2 >= 1000000007. (2 <= 38)
 |  144004141^38 mod 1000000007 = 144004141^{2*19+0} mod 1000000007 = 144004141^{2*19} * 144004141^0 mod 1000000007
 |  = (144004141^2 mod 1000000007)^19 * 144004141^0 mod 1000000007
 |  = 479987537^19 * 144004141^0 mod 1000000007
 |  = 17^19 * 131^19 * 215531^19 mod 1000000007
 |  Then we divide the expression in 3 parts.
 | |  17^8 >= 1000000007. (8 <= 19)
 | |  17^19 mod 1000000007 = 17^{8*2+3} mod 1000000007 = 17^{8*2} * 17^3 mod 1000000007
 | |  = (17^8 mod 1000000007)^2 * 17^3 mod 1000000007
 | |  = 975757399^2 * 17^3 mod 1000000007
 | |  = 17^3 * 59^2 * 16538261^2 mod 1000000007
 | |  Then we divide the expression in 3 parts.
 | | |  16538261^2 >= 1000000007. (2 <= 2)
 | | |  16538261^2 mod 1000000007 = 16538261^{2*1+0} mod 1000000007 = 16538261^{2*1} * 16538261^0 mod 1000000007
 | | |  = (16538261^2 mod 1000000007)^1 * 16538261^0 mod 1000000007
 | | |  = 74989523^1 * 16538261^0 mod 1000000007
 | | |  = 7 * 19 * 563831 mod 1000000007
 | | |  Then we divide the expression in 3 parts.
 | | |  7 * 19 * 563831 mod 1000000007 = 74989523
 | | |

 | |  4913 * 3481 * 74989523 mod 1000000007 = 286765645
 | |

 | |  131^5 >= 1000000007. (5 <= 19)
 | |  131^19 mod 1000000007 = 131^{5*3+4} mod 1000000007 = 131^{5*3} * 131^4 mod 1000000007
 | |  = (131^5 mod 1000000007)^3 * 131^4 mod 1000000007
 | |  = 579489385^3 * 131^4 mod 1000000007
 | |  = 5^3 * 131^4 * 115897877^3 mod 1000000007
 | |  Then we divide the expression in 3 parts.
 | | |  115897877^2 >= 1000000007. (2 <= 3)
 | | |  115897877^3 mod 1000000007 = 115897877^{2*1+1} mod 1000000007 = 115897877^{2*1} * 115897877^1 mod 1000000007
 | | |  = (115897877^2 mod 1000000007)^1 * 115897877^1 mod 1000000007
 | | |  = 799080910^1 * 115897877^1 mod 1000000007
 | | |  = 2 * 5 * 19 * 4205689 * 115897877 mod 1000000007
 | | |  Then we divide the expression in 5 parts.
 | | |  2 * 5 * 19 * 4205689 * 115897877 mod 1000000007 = 371945610
 | | |

 | |  125 * 294499921 * 371945610 mod 1000000007 = 316393257
 | |

 | |  215531^2 >= 1000000007. (2 <= 19)
 | |  215531^19 mod 1000000007 = 215531^{2*9+1} mod 1000000007 = 215531^{2*9} * 215531^1 mod 1000000007
 | |  = (215531^2 mod 1000000007)^9 * 215531^1 mod 1000000007
 | |  = 453611639^9 * 215531^1 mod 1000000007
 | |  = 13^9 * 215531 * 34893203^9 mod 1000000007
 | |  Then we divide the expression in 3 parts.
 | | |  13^9 >= 1000000007. (9 <= 9)
 | | |  13^9 mod 1000000007 = 13^{9*1+0} mod 1000000007 = 13^{9*1} * 13^0 mod 1000000007
 | | |  = (13^9 mod 1000000007)^1 * 13^0 mod 1000000007
 | | |  = 604499303^1 * 13^0 mod 1000000007
 | | |  = 37 * 89 * 183571 mod 1000000007
 | | |  Then we divide the expression in 3 parts.
 | | |  37 * 89 * 183571 mod 1000000007 = 604499303
 | | |

 | | |  34893203^2 >= 1000000007. (2 <= 9)
 | | |  34893203^9 mod 1000000007 = 34893203^{2*4+1} mod 1000000007 = 34893203^{2*4} * 34893203^1 mod 1000000007
 | | |  = (34893203^2 mod 1000000007)^4 * 34893203^1 mod 1000000007
 | | |  = 607076464^4 * 34893203^1 mod 1000000007
 | | |  = 2^16 * 37^4 * 733^4 * 1399^4 * 34893203 mod 1000000007
 | | |  Then we divide the expression in 5 parts.
 | | | |  733^4 >= 1000000007. (4 <= 4)
 | | | |  733^4 mod 1000000007 = 733^{4*1+0} mod 1000000007 = 733^{4*1} * 733^0 mod 1000000007
 | | | |  = (733^4 mod 1000000007)^1 * 733^0 mod 1000000007
 | | | |  = 679467505^1 * 733^0 mod 1000000007
 | | | |  = 5 * 8543 * 15907 mod 1000000007
 | | | |  Then we divide the expression in 3 parts.
 | | | |  5 * 8543 * 15907 mod 1000000007 = 679467505
 | | | |

 | | | |  1399^3 >= 1000000007. (3 <= 4)
 | | | |  1399^4 mod 1000000007 = 1399^{3*1+1} mod 1000000007 = 1399^{3*1} * 1399^1 mod 1000000007
 | | | |  = (1399^3 mod 1000000007)^1 * 1399^1 mod 1000000007
 | | | |  = 738124185^1 * 1399^1 mod 1000000007
 | | | |  = 3 * 5 * 1399 * 49208279 mod 1000000007
 | | | |  Then we divide the expression in 4 parts.
 | | | |  3 * 5 * 1399 * 49208279 mod 1000000007 = 635727591
 | | | |

 | | |  65536 * 1874161 * 679467505 * 635727591 * 34893203 mod 1000000007 = 435840952
 | | |

 | |  604499303 * 215531 * 435840952 mod 1000000007 = 400832125
 | |

 |  286765645 * 316393257 * 400832125 mod 1000000007 = 685365020
 |

  877905026 * 2211169 * 685365020 mod 1000000007 = 779596752


 1048576 * 390957170 * 580316551 * 779596752 mod 1000000007 = 676888799

Simple Looping

Of course, another simple and navie solution is to loop power times.

public static int FindModularExponentiationLoop(int b, int p, int m)
{
	int result = 1;
	while (p-- >= 1)
		result = (int) ((long)result * b % m);
	return result;
}

The following is the performance evaluation by BenchmarkDotNet. Unfortunately my "smart" method is terribly slow.

Method	Mean	Error	StdDev
Smart	299.70 us	0.446 us	0.396 us
Loop	47.46 us	0.081 us	0.075 us

参考资料

Will Calderwood. Am I going crazy or is Math.Pow broken?. . 2011-03-23 [2024-05-20].↑
SunsetQuest. How do you do *integer* exponentiation in C#?. . 2014-02-13 [2024-05-20].↑

《Code for 大数取模 modulo arithmetic》有一个想法

大数取模 modulo arithmetic – gqqnbig的专栏说道：

2024年5月19日 21:11

[…] Up to now, we have to use a computer. See Code for 大数取模 modulo arithmetic. […]

评论已关闭。