code: def qs(N, prime_bound, M, ERROR_TERM=25, seed=1234):
      """Performs factorization using Self-Initializing Quadratic Sieve.
      In SIQS, let N be a number to be factored, and this N should not be a
      perfect power. If we find two integers such that ``X**2 = Y**2 modN`` and
      ``X != +-Y modN``, then `gcd(X + Y, N)` will reveal a proper factor of N.
      In order to find these integers X and Y we try to find relations of form
      t**2 = u modN where u is a product of small primes. If we have enough of
      these relations then we can form ``(t1*t2...ti)**2 = u1*u2...ui modN`` such that
      the right hand side is a square, thus we found a relation of ``X**2 = Y**2 modN``.
  
      Here, several optimizations are done like using multiple polynomials for
      sieving, fast changing between polynomials and using partial relations.
      The use of partial relations can speeds up the factoring by 2 times.
  
      Parameters
      ==========
  
      N : Number to be Factored
      prime_bound : upper bound for primes in the factor base
      M : Sieve Interval
      ERROR_TERM : Error term for checking smoothness
      threshold : Extra smooth relations for factorization
      seed : generate pseudo prime numbers
  
      Examples
      ========
  
      >>> from sympy.ntheory import qs
      >>> qs(25645121643901801, 2000, 10000)
      {5394769, 4753701529}
      >>> qs(9804659461513846513, 2000, 10000)
      {4641991, 2112166839943}
  
      References
      ==========
  
      .. [1] https://pdfs.semanticscholar.org/5c52/8a975c1405bd35c65993abf5a4edb667c1db.pdf
      .. [2] https://www.rieselprime.de/ziki/Self-initializing_quadratic_sieve
      """
      ERROR_TERM *= 2**10
      idx_1000, idx_5000, factor_base = _generate_factor_base(
          prime_bound, N
      )
      smooth_relations = []
      ith_poly = 0
      partial_relations = {}
      proper_factor = set()
      threshold = 5 * len(factor_base) // 100
      while True:
          if ith_poly == 0:
              ith_sieve_poly, B_array = _initialize_first_polynomial(
                  N, M, factor_base, idx_1000, idx_5000
              )
          else:
              ith_sieve_poly = _initialize_ith_poly(
                  N, factor_base, ith_poly, ith_sieve_poly, B_array
              )
          ith_poly += 1
          if ith_poly >= 2 ** (
              len(B_array) - 1
          ):  # time to start with a new sieve polynomial
              ith_poly = 0
          sieve_array = _gen_sieve_array(M, factor_base)
          s_rel, p_f = _trial_division_stage(
              N,
              M,
              factor_base,
              sieve_array,
              ith_sieve_poly,
              partial_relations,
              ERROR_TERM,
          )
          smooth_relations += s_rel
          proper_factor |= p_f
          if len(smooth_relations) >= len(factor_base) + threshold:
              break
      matrix = _build_matrix(smooth_relations)
      dependent_row, mark, gauss_matrix = _gauss_mod_2(matrix)
      N_copy = N
      for index in range(len(dependent_row)):
          factor = _find_factor(
              dependent_row,
              mark,
              gauss_matrix,
              index,
              smooth_relations,
              N,
          )
          if factor > 1 and factor < N:
              proper_factor.add(factor)
              while N_copy % factor == 0:
                  N_copy //= factor
              if isprime(N_copy):
                  proper_factor.add(N_copy)
                  break
              if N_copy == 1:
                  break
      return proper_factor
  
imports: from sympy.core.random import _randint
  from sympy.external.gmpy import gcd, invert, sqrt
  from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power
  from sympy.ntheory import isprime
  from math import log, sqrt
label_desc: Use image node_id 10 for calling a global function with example usage: qs(N, prime_bound, M, ERROR_TERM, seed) and returns: proper_factor