colt_browning

negation & addition of Gaussian integers written in base i-1

Oct 12th, 2016
278
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. def neg(n):
  2.     '''Negate a Khmelnik-Penney-encoded Gaussian integer.
  3.  
  4. n is interpreted as a binary code with base i-1,
  5. where i is the imaginary unit.
  6. Algorithm from Solomon I. Khmelnik,
  7. "Specialized Digital Computer for Operations with Complex Numbers"
  8. (in Russian), Questions of Radio Electronics, volume 12, number 2, pp. 60-82, 1964,
  9. http://lib.izdatelstwo.com/Papers2/s4.djvu
  10. See table 4.
  11.  
  12. See also oeis.org/A066321, references therein, and the related SeqFan threads:
  13. http://list.seqfan.eu/pipermail/seqfan/2016-August/016646.html
  14. http://list.seqfan.eu/pipermail/seqfan/2016-August/016650.html
  15. http://list.seqfan.eu/pipermail/seqfan/2016-September/016716.html
  16. http://list.seqfan.eu/pipermail/seqfan/2016-October/016843.html'''
  17.     sigma, shift = 0, 0
  18.     while n >> shift:
  19.         bits = n&(1<<shift)
  20.         shift += 1
  21.         sigma |= bits
  22.         if bits:
  23.             bits = n&(1<<shift)
  24.             shift += 1
  25.             sigma |= bits
  26.             if bits:
  27.                 bits = ~n&(1<<shift)
  28.                 shift += 1
  29.             else:
  30.                 bits = ~n&(7<<shift)
  31.                 shift += 3
  32.             sigma |= bits
  33.     return sigma
  34.  
  35. def neg2(n):
  36.     '''Negate a Khmelnik-Penney-encoded Gaussian integer.
  37.  
  38. A shorter, more cryptic and slightly slower version of neg(n).'''
  39.     n0, sigma, shift = n, 0, 0
  40.     while n:
  41.         if n&1:
  42.             sigma |= (~n & 31>>(n&2) ^ 3) << shift
  43.         shift += [1, 5, 1, 3][n%4]
  44.         n = n0 >> shift
  45.     return sigma
  46.  
  47. def add(a, b):
  48.     '''Add up two Khmelnik-Penney-encoded Gaussian integers.
  49.  
  50. a and b are interpreted as binary codes with base i-1.
  51. Algorithm from Solomon I. Khmelnik,
  52. "Specialized Digital Computer for Operations with Complex Numbers"
  53. (in Russian), Questions of Radio Electronics, volume 12, number 2, pp. 60-82, 1964,
  54. http://lib.izdatelstwo.com/Papers2/s4.djvu
  55. See table 7.'''
  56.     sigma, shift, sk = 0, 0, 0o4444
  57.     # Every octal digit of sk corresponds to one of the "retarded carries" P.
  58.     while a or b or sk not in [0o4444, 0o4322]:
  59.         sk = ((sk>>3) + [0o5000, 0o3670, 0o4000, 0o2670][(sk>>1)%4]
  60.               + (a&1) + (b&1))
  61.         sigma |= (sk&1) << shift
  62.         shift += 1
  63.         a >>= 1
  64.         b >>= 1
  65.     return sigma
  66.  
  67. def from_khmelnik(n, cx = False):
  68.     '''Obtain a Khmelnik-Penney-encoded Gaussian integer.
  69.  
  70. Returns a pair of integers, i. e., real and imaginary
  71. parts of the result. Alternatively, if the second argument
  72. is set to True, returns a complex number.'''
  73.     gre, gim, bre, bim = 0, 0, 1, 0
  74.     while n:
  75.         if n&1:
  76.             gre += bre
  77.             gim += bim
  78.         n >>= 1
  79.         bre, bim = -bre-bim, bre-bim
  80.     return gre+1j*gim if cx else (gre, gim)
  81.  
  82. def to_khmelnik(re, im = None):
  83.     '''Obtain Khmelnik-Penney encoding of a Gaussian integer.
  84.  
  85. Accepts either a complex number or a pair of integers
  86. (real and imaginary parts).'''
  87.     if im is None:
  88.         if isinstance(re, int):
  89.             im = 0
  90.         else:
  91.             re, im = round(re.real), round(re.imag)
  92.     n = 0 # to_khmelnik(0)
  93.     k = 1 # to_khmelnik(1)
  94.     if re < 0:
  95.         k = neg(k)
  96.         re = -re
  97.     while re:
  98.         if re&1:
  99.             n = add(n, k)
  100.         re >>= 1
  101.         k = add(k, k)
  102.     k = 3 # to_khmelnik(1j)
  103.     if im < 0:
  104.         k = neg(k)
  105.         im = -im
  106.     while im:
  107.         if im&1:
  108.             n = add(n, k)
  109.         im >>= 1
  110.         k = add(k, k)
  111.     return n
RAW Paste Data