PNlibSrc = md`\`\`\`python """\ PolyNum class ============ (Polynomial Number class) Floating point arithmetic for "numbers" with digit like real number or items from any arbitrary field, allow solving differential or difference equations using Mikusinski's or Bellert's approach. https://pei.prz.edu.pl/%7Ekubaszek/index_en.html We assume PN in form: (~m_0 ~, m_1 ~ m_2 ~ ...~) · (~1~0~)^c where: m - PN mantissa c - PN exponent version: 0.30, 2019-04-16, Brython todo: more convenient PN.exp() if self.exponent > 0 """ #2018-05-24 x is scalar == isinstance(x,(int,float)) or not hasattr(x, '__len__') #__debugger__ PolyNumConf_max_N = 64 #or 128 or 100 ..., 32 is mimnimum for doctest - PN significant digits number PolyNumConf_sep = '~' # in PN str() and repr() #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@# from browser import window epsilonPNdig = window.Number.EPSILON #sys.float_info.epsilon #=2**-52 = 2.220446049250313e-16 flt = float # def flt(digStr): return(float(digStr)) def strF(dig, chop_=True, signifi_=9): if dig and chop_: dig = chop(dig) return dig.__format__('1.'+str(signifi_)) # try: return dig.__format__('1.'+str(signifi_)) #'1.9' # except (ValueError, TypeError): return str(dig) # #without formatting (i.e. for int, TypeError - for Fraction) # universal ''reprF'', ''zeroPNdig'', ''onePNdig'' and ''PNdig_isclose'' def reprF(dig): """ repr() of PN digit """ # return repr(dig) return strF(dig, chop_=False, signifi_=15) zeroPNdig = flt('0') # flt('0') #zero of type of PN digits onePNdig = flt('1') def chop(d, tol=epsilonPNdig*1024*1024, zero=zeroPNdig):#tol: 6e-33 by eps=6e-39 """ Converts x close to zero to exact zero """ if (abs(d) > tol): return d else: return zero def PNdig_isclose(a, b, rel_tol=epsilonPNdig*128, abs_tol=epsilonPNdig*128): """ Like ''math.isclose'' or ''mpmath.almosteq'' - for an arbitrary types a, b Return : True if a - b == 0 or |a - b| <= abs_tol else if (|a-b| and rel_tol) are compatible type: |a-b| <= rel_tol*max(|a|,|b|) else: raise NotImplementedError """ if not a and not b: return True diff = abs(a - b) if not diff: return True if diff <= abs_tol: return True absMax, abs2 = abs(a), abs(b) if absMax < abs2: absMax = abs2 if isinstance(diff,type(rel_tol)): return diff <= rel_tol * absMax raise NotImplementedError #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@# floor = window.Math.floor sqrt = window.Math.sqrt exp = window.Math.exp log = window.Math.log pi = window.Math.PI #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@# class PolyNum(object): """\ A Polynomial Number class ========================= Floating point arithmetic for "numbers" with digit like real number or items from any arbitrary field, allow solving differential or difference equations using Mikusinski's or Bellert's approach. https://pei.prz.edu.pl/%7Ekubaszek/index_en.html We assume **Polynomial Number** (**PN**) in form: (~m_0 ~, m_1 ~ m_2 ~ ...~) · (~1~0~)^c where: m - PN mantissa c - PN exponent Attributes: ---------- mantissa: the PN "digits", in decreasing (negative) powers, starting from power 0. List of _max_N elements, m[0] * p**(0) + m[1] * p**(-1) + ... +m[_max_N-1] * p**(-(_max_N-1)) If initial length is shorten then _max_N, mantissa is filled with zeros of type as first element. exponent: (int) _strPN_cut: (int) (not defined by default) could be added to istance, to show first _strPN_cut digits in str() and repr(). Parameters ---------- mantissa_or_pN_or_str : array_like type: the Polynomial Number's coefficients, in decreasing (negative) powers, starting from power 0. (see class attr. self._max_N) PN type: PN to copy on __init__, adding exponentAdd to exponent string type: '*(~1~0~)**()' where is str(integer), is '(~~,~~...~)' ,... - real numbers, '~,' (or '~,~') denotes "radix point" (see class attr. self._sep = '~') 'const:(~2~,-4~4~-4~4~...~)' p_trap(h)==(2/h)*(~1~-1~)/(~1~1~)==(1/h)*(~2~,-4~4~-4~4~...~) 'const:(~1~,2~2~2~2~2~2~...~)' flat samples ready to multiply by (a/2): exponentAdd: value added to 0 or original PN exponent (if exists) Const. fixed as class attr. --------------------------- _max_N = PolyNumConf_max_N """ #@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@# __hash__ = None _max_N = PolyNumConf_max_N #num. of mantissa items _sep = PolyNumConf_sep #'~' # in str() and repr() @property def mantissa(self): """ A copy of the coefficients""" return self._mantissa[:] #or deepcopy ? @property def exponent(self): """ The exponent of the Polynomial Number """ if hasattr(self, '_exponent'): # ? return self._exponent else: return 0 def __init__(self, mantissa_or_pN_or_str=None, exponentAdd=0): """\ After __init__ first PN digit (self.mantissa[0]) is non-zero or (in zero PN case) mantissa is filled with zeroes. """ if isinstance(mantissa_or_pN_or_str, PolyNum): self._initMant(mantissa_or_pN_or_str.mantissa) self._exponent = exponentAdd if hasattr(self,'exponent'): self._exponent = mantissa_or_pN_or_str.exponent + exponentAdd return if mantissa_or_pN_or_str is None: mantissa_or_pN_or_str = [] # elif isinstance(mantissa_or_pN_or_str, string_types): # six.string_types elif isinstance(mantissa_or_pN_or_str, str): #__future__ unicode_literals mantissa_or_pN_or_str, exponent1 = \ getMantissaExponent_fromStr(mantissa_or_pN_or_str, self._sep, self._max_N) exponentAdd += exponent1 else: #mantissa_or_pN_or_str should be array_like #if isinstance(mantissa_or_pN_or_str, numbers.Number): if not isinstance(mantissa_or_pN_or_str,(int,float)) or \ hasattr(mantissa_or_pN_or_str, '__len__'): mantissa_or_pN_or_str = list(mantissa_or_pN_or_str) else: #scalar mantissa_or_pN_or_str = [mantissa_or_pN_or_str] self._initMant(mantissa_or_pN_or_str) self._exponent = exponentAdd self._normalize0() #case: PolyNum('(~0~,2.~-0.3~)'), 100 ## print('init>>>',self) def _initMant(self, mant=None): """ default init: [0.0, 0.0, ...] zeros of type zeroPNdig """ if isinstance(mant, PolyNum): self._mantissa = mant.mantissa return #if not mant: #numpy ValueError: The truth value of an array with more than one element is ambiguous. if (mant is None) or not len(mant): mant = [zeroPNdig] #type is determined in digitPN, rather not [0.0] self._mantissa = list(mant[:self._max_N]) self._fillZerosProc() def _fillZerosProc(self): assert self._mantissa, "At least one element must exist, to add zeroes of the same type!" if len(self._mantissa) < self._max_N: zero = 0 * self._mantissa[0] self._mantissa += [zero for __ in range(self._max_N-len(self._mantissa))] def _shrMantProc(self, r): zero = 0 * self._mantissa[0] if r > self._max_N: r = self._max_N zeros = [zero for __ in range(r)] self._mantissa = zeros + self._mantissa[:(self._max_N - r)] def _normalize0(self): #return None """First digit non-zero or zero PN - to use in init""" if self._mantissa[0]: #nothing to do return self first = 0 m0 = self._mantissa[0] for d in self._mantissa: if d: #d!= 0 break else: first += 1 self._exponent -= 1 if first == len(self._mantissa): self._mantissa = [m0] #preserve first 0.0 - to get type of mant[0] else: self._mantissa = self._mantissa[first:] self._fillZerosProc() # fill right zeros #----------------------------------------------------------- def _normalize(self): """ return PolyNum copy with first digit non-zero or zero PN if actions with mantissa/exponent was performed. If first digit is nonzero return self reference (not a copy). - To use after add, sub, div """ if self._mantissa[0]: #nothing to do return self y = PolyNum(self) #copy of mantissa, because it will be shifted y._normalize0() return y def __str__(self): ''' Examples -------- (~1.2~,2.5~-0.3~)*(~1~0~)**(-2) (~1.2~,2.5~-0.3~)*(~1~0~)**(2) (~1.2~,2.5~-0.3~) (~1.2~) ''' return self._strPN(strF) def __repr__(self): return "PolyNum('"+self._strPN(reprF)+"')" def _strPN(self, digToStr): """ digToStr = strF """ cutedStr = hasattr(self, '_strPN_cut') and self._strPN_cut < self._max_N ma = self.mantissa ex = self.exponent # r = '({0}{1}{0}'.format(self._sep, ma[0]) # (~{}~ try: r1 = digToStr(ma[0]) except ValueError: #try without formatting (i.e. for int) r1 = str(ma[0]) r = '('+self._sep + r1 + self._sep #~ if ma[0]: #non-zero PN ma = ma[1:] #ma = trim_zeros(ma, trim='b') #'b' to trim from back last = len(ma) #for d in ma[::-1]: #transcrypt doesn''t like ''[::-1]'' for d in reversed(ma): if d: break else: last = last - 1 ma = ma[:last] if cutedStr and (len(ma) < self._strPN_cut-1): cutedStr = False if cutedStr: ma = ma[:self._strPN_cut-1] if len(ma): r += ',' for m in ma: r += digToStr(m) + self._sep #~ if cutedStr: r += '...~' r += ')' if ex and self.__nonzero__(): r += '*({0}1{0}0{0})**({1})'.format(self._sep, ex) #*(~1~0~)**({}) return r def chop(self, tol=epsilonPNdig*1024*1024): #1024*1024 = 1e6 """ Converts PN digits close to zero to exact zeros. Normalize at the end. todo: tol can grow with digit index """ zero = self.mantissa[0] * 0 ma = [chop( d, tol, zero ) for d in self.mantissa] x = PolyNum(ma, self.exponent) return x._normalize0() def __format__(self, format_spec): s = self.__str__() return s.__format__(format_spec) def __nonzero__(self): return self._mantissa[0] def __len__(self): return self._max_N def __call__(self, val): """ """ if not self.__nonzero__(): return val * 0 # zero of val type y = mantPN_val(self.mantissa, val) res = 1 if self.exponent: res = val #val **self.exponent for _ in range(abs(self.exponent)-1): res *= val if self.exponent < 0: return y / res else: return y * res def __neg__(self): return PolyNum([-x for x in self._mantissa], self.exponent) def __pos__(self): return self def __mul__(self, other): """ """ if isinstance(other, PolyNum): yMant = mantPN_mul(self.mantissa, other.mantissa, self._max_N) elif isinstance(other,(int,float)) or \ not hasattr(other, '__len__'): #''other'' is scalar yMant = [x * other for x in self.mantissa] else: raise ValueError(str(other)+" <- invalid operand for PN __mul__") expoOther = 0 if hasattr(other, 'exponent'): expoOther = other.exponent return PolyNum(yMant, self.exponent + expoOther) def __rmul__(self, other): # case: other * self if isinstance(other, PolyNum): yMant = mantPN_mul(other.mantissa, self.mantissa, self._max_N) elif isinstance(other,(int,float)) or \ not hasattr(other, '__len__'): #''other'' is scalar yMant = [other * x for x in self.mantissa] else: raise ValueError(str(other)+" <- invalid operand for PN __rmul__") expoOther = 0 #f.ex. int, real, if hasattr(other, 'exponent'): expoOther = other.exponent return PolyNum(yMant, self.exponent + expoOther) def __truediv__(self, other): #self / other """ """ if isinstance(other, PolyNum): yMant = mantPN_mul( self.mantissa, mantPN_inv(other.mantissa, self._max_N), self._max_N) elif isinstance(other,(int,float)) or \ not hasattr(other, '__len__'): #''other'' is scalar yMant = [x / other for x in self.mantissa] else: raise ValueError(str(other)+" <- invalid operand for PN __div__") expoOther = 0 if hasattr(other, 'exponent'): expoOther = other.exponent return PolyNum(yMant, self.exponent - expoOther) #? __div__ = __truediv__ def __rtruediv__(self, other): # case: other / self inv = mantPN_inv(self.mantissa, self._max_N) if isinstance(other, PolyNum): yMant = mantPN_mul(other.mantissa, inv, self._max_N) else: if other == 1 or other == onePNdig: yMant = inv elif isinstance(other,(int,float)) or \ not hasattr(other, '__len__'): #''other'' is scalar yMant = [other * x for x in inv] else: raise ValueError(str(other)+" <- invalid operand for PN __rdiv__") expoOther = 0 if hasattr(other, 'exponent'): expoOther = other.exponent return PolyNum(yMant, expoOther - self.exponent) #? __rdiv__ = __rtruediv__ def __idiv__(self, other): # need for transcript return self.__truediv__(other) def __rdiv__(self, other): # need for transcript return self.__rtruediv__(other) def __div__(self, other): #? return self.__truediv__(other) def __add__(self, other): """ """ if ( isinstance(other,(int,float)) or \ not hasattr(other, '__len__') ) and not other: return PolyNum(self) #PolyNum() - copy of mantissa otherPN = PolyNum(other) if not self.__nonzero__(): return PolyNum(otherPN) if otherPN.exponent <= self.exponent: a1, a2 = self, otherPN else: a1, a2 = otherPN, self if a2.exponent < a1.exponent - self._max_N: return PolyNum(a1) if a2.exponent == a1.exponent: #it will be a common case - do not shrMant() # yMant = [a + b for (a,b) in zip_longest(a1.mantissa, a2.mantissa, fillvalue=0)]) yMant = [a + b for (a,b) in zip(a1.mantissa, a2.mantissa)] else: a2 = PolyNum(a2) #copy of a2 mantissa, which will be shifted a2._shrMantProc(a1.exponent - a2.exponent) #because of leading zero a2 is treating as (~0~) PolynNum in debugger, #but i.e. a2.mantissa == [0.0, 0.0, 1.0, 2, 3, 0.0, ... yMant = [a + b for (a,b) in zip(a1.mantissa, a2.mantissa)] return PolyNum(yMant, a1.exponent)._normalize() def __radd__(self, other): return self.__add__(PolyNum(other)) def __sub__(self, other): ''' ''' return self.__add__(-PolyNum(other)) def __rsub__(self, other): # case: other - self selfNeg = -self return selfNeg.__radd__(other) def __pow__(self, a): """ self **a ======== where exists self[0] **a (a is rather float, fract or int) """ if not a or not self.__nonzero__(): # results of x[0] **0 and 0 **a (rather =1) determined by mantissa[0] **a return PolyNum([self._mantissa[0] **a]) # (~1~,0~0~...~) #if (self.exponent != 0 and not isinstance(a, numbers.Integral)): if (self.exponent != 0 and not isinstance(a, int)): raise ValueError("Power only to int, real, rational by exponent == 0 or to int by exponent != 0") #if self.exponent == 0 power to int, real, rational allowed yMant = mantPN_power_real(self.mantissa, a, self._max_N) expo = 0 if self.exponent != 0: expo = self.exponent * a #assert isinstance(expo, numbers.Integral) assert isinstance(expo, int) return PolyNum(yMant, expo) def _MantPN_isclose(self, other, rel_tol=epsilonPNdig*128, abs_tol=epsilonPNdig*128): """ Return : True if True for all pair digits a,b from self and others: if a - b == 0 or |a - b| <= abs_tol: True else: if (|a-b| and rel_tol) are compatybile type: |a-b| <= rel_tol*max(|a|,|b|) else: raise NotImplementedError """ if not isinstance(other, PolyNum): return NotImplemented for a, b in zip(self.mantissa, other.mantissa): if not PNdig_isclose(a, b, rel_tol, abs_tol): # a != b: return False return True def isclose(self, other, rel_tol=epsilonPNdig*128, abs_tol=epsilonPNdig*128): """ Return : True if True for all pair digits a,b from self and others: if a - b == 0 or |a - b| <= abs_tol: True else: if (|a-b| and rel_tol) are compatybile type: |a-b| <= rel_tol*max(|a|,|b|) else: raise NotImplementedError """ if (not self) and (not other): return True if isinstance(other, PolyNum): otherPN = other else: otherPN = PolyNum(other) if otherPN.exponent <= self.exponent: a1, a2 = self, otherPN else: a1, a2 = otherPN, self if a2.exponent != a1.exponent: a2 = PolyNum(a2) #copy of a2 mantissa, which will be shifted a2._shrMantProc(a1.exponent - a2.exponent) #because of leading zero a2 is treating as (~0~) PolynNum in debugger, #but i.e. a2.mantissa == [0.0, 0.0, 1.0, 2, 3, 0.0, ... return a1._MantPN_isclose(other, rel_tol, abs_tol) def __eq__(self, other): return self.isclose(other) def __ne__(self, other): if not isinstance(other, PolyNum): return NotImplemented return not self.__eq__(other) def __abs__(self): """ | . | see a concept of partially ordered Mikusinski space """ return PolyNum([abs(x) for x in self.mantissa], self.exponent) def isnonnegative(self): """ self >= 0 see a concept of partially ordered Mikusinski space """ if not self.__nonzero__(): return True for x in self.mantissa: if not x >= 0: return False return True def __le__(self, other): """ self <= other see a concept of partially ordered Mikusinski space some PNs are not comparable, i.e. ~(self <= other) and ~(self >= other) """ other = PolyNum(other) return ( other - self ).isnonnegative() def __ge__(self, other): """ self >= other see a concept of partially ordered Mikusinski space some PNs are not comparable, i.e. ~(self >= other) and ~(self <= other) """ return ( self - other ).isnonnegative() ######## def sqrt(self): """ """ if self._exponent and (self.exponent % 2): #odd raise ValueError("Does not support sqrt() if exponent is odd: {}.format()self.exponent") return PolyNum(mantPN_sqrt(self.mantissa, self._max_N), self._exponent // 2) def exp(self): """ exp() ======== >>> y3 = PolyNum('(~1.~2.~0.)').exp() Traceback (most recent call last): ... ValueError: Todo(?) - case if exponent (== 2) is positive. """ if self._exponent > 0: raise ValueError("Todo(?) - case if exponent (== {}) is positive.".format(self.exponent)) if self._exponent == 0: #it will be common case - do not shrMant() return PolyNum(mantPN_exp(self.mantissa, self._max_N)) else: a2 = PolyNum(self) #copy of a2 mantissa, which will be shifted a2._shrMantProc(-a2._exponent) # because of leading zero a2 is treating as (~0~), but i.e. # a2.mantissa [0.0, 0.0, 1.0, 2, 3, 0.0, ... return PolyNum(mantPN_exp(a2.mantissa, self._max_N)) #### def expZ(self, pZ, T0, h): """ returns exp() with shifted samples exp(arg+K·h·pZ) * (~1~0~)^(-K) where K = floor( T0 / h ) T0 - expected pure delay time h - sampling interval pZ - numerical Heaviside operator, e.g. 2/h * (~1~-1~)/(~1~1~) """ K = int(floor( T0 / h )) if not K: return self.exp() else: x = self + K*h*pZ return PolyNum(x.exp(), -K) #### def ln(self): if self._exponent: raise ValueError("Does not support ln() if exponent is nonzero: {}.format()self.exponent") return PolyNum(mantPN_ln(self.mantissa, self._max_N)) def __getitem__(self, index): """ index < 0 is not allowed index == i + (-exponent) exponent > 0 is not allowed useful for Z-transform, like mantissa normalized to exponent==0, where exponent <= 0; exponent > 0 no allowed used in plot """ if not self.exponent: #the most common case return self.mantissa[index] if type(index) is slice: if self._exponent > 0: raise ValueError("if exponent(={}) is positive, PN slicing is not allowed".format(self._exponent)) y = ([0 * self._mantissa[0]] * (-self._exponent)) \ + self.mantissa[:self._exponent] # [:-(-self._exponent)] return y[index] else: # int if not 0 <= index < len(self._mantissa): raise IndexError() i = index + self._exponent #if i < 0 : raise IndexError() # ^ ValueError: cannot copy sequence with size 32 to array axis with dimension 0 #if i >= self._max_N: raise IndexError() # ^ ValueError: cannot copy sequence with size 32 to array axis with dimension 33 # ? if exp == 2 # ValueError: cannot copy sequence with size 32 to array axis with dimension 30 if i < 0: return 0 * self._mantissa[0] # zero of type m[0] else: return self.mantissa[i] # def __setitem__(self, index, val): # self._mantissa[index] = val # return def __iter__(self): """ if __iter__ is not defined, then __getitem__ is used for list() """ if not self.exponent: #the most common case return iter(self.mantissa) if self._exponent > 0: raise ValueError("if exponent(={}) is positive, PN __iter__ is not allowed".format(self._exponent)) zero = 0 * self._mantissa[0] y = ([zero for __ in range(-self._exponent)]) + self.mantissa[:self._exponent] # [:-(-self._exponent)] return iter(y) def asList(self): ''' List normalized to exponent=0 - only if exponent <= 0 Can be longer than _max_N ''' if not self.__nonzero__: mantissaE0 = [self._mantissa[0] for __ in range(self._max_N)] #[0.,0.,0.,...] else: ex = -self.exponent if ex == 0: mantissaE0 = self.mantissa elif ex < 0: raise ValueError( 'PN exponent = {}. Can not convert to array if PN exponent > 0'.format(-ex)) else: #(-self.exponent) > 0 zero = 0 * self._mantissa[0] mantissaE0 = [zero for __ in range(ex)] + self._mantissa return mantissaE0 def asPointsList(self, h, t_offs=0): ''' [[x,y],...]; x=i·h; y.asList() - for plot ''' return [[i * h + t_offs,y] for i, y in enumerate(self.asList())] #import numpy as np #return np.asarray(mantissaE0) # def __array__(self): return np.asarray(self.asList) #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& def invTr1LaplPN (self, t): """ inp: t :float otpt: y, err :float Examples: >>> Yerr = PolyNum('(~1~2~)') >>> Yerr.invTr1LaplPN(1) Traceback (most recent call last): ... ValueError: PN exponent = 1 > 0 - inverese Laplce transform does not exist. """ stop = False zero = self._mantissa[0] * 0 one = zero + 1 #1 of type of zero if t < 0: return zero, zero max_a = zero ex = self._exponent if ex >= 0: raise ValueError( 'PN exponent = {} > 0 - inverese Laplce transform does not exist.'.format(ex)) if t == 0: #return self._mantissa[-ex-1], zero #rather set t=0+ = eps/int(1e24) t = epsilonPNdig / 1000000000000000000000000 wk = one for k in range(1, -ex-1+1): # (~0~,0~0~0~0 wk = wk * t/k out= self._mantissa[0] * wk abs_a = [abs(out),zero,zero] #ring-array of 3 last abs_ak[0..2] dk_maxLocal = -1 k = 0 while 1: wk = wk * t / (k-ex) ak = self._mantissa[k+1] * wk abs_ak = abs(ak) if (abs_ak > max_a): max_a = abs_ak abs_a[(k+1) % 3] = abs_ak if (abs_a[(k+3-1) % 3] < abs_a[k % 3]) \ and (abs_a[k % 3] >= abs_a[(k+1) % 3]): # local max is here if dk_maxLocal < 0: #start maxLocal = abs_a[k % 3] k_maxLocal = k dk_maxLocal = 0 else: q_maxLocal = abs_a[k % 3]/maxLocal maxLocal = abs_a[k % 3] dk_maxLocal = k-k_maxLocal k_maxLocal = k stop = ( (dk_maxLocal > 0) and (q_maxLocal < 1) and (maxLocal < abs(out)*epsilonPNdig) ) \ or \ ( (k > k_maxLocal + self._max_N // 4) and (abs_ak < abs(out)*epsilonPNdig) ) #last line - not oscillatory convergent out=out + ak k += 1 if (k >= self._max_N-1) or stop: break if ( (k <= k_maxLocal+self._max_N // 4) and (dk_maxLocal > 0) and (q_maxLocal < 1) ): err = maxLocal * q_maxLocal **((k-k_maxLocal)/dk_maxLocal ) else: err = abs_ak err = err + (max_a+out) * epsilonPNdig; return out, err #end def invTr1LaplPN(self, t) def invTr05exp_b0_LaplPN (self, t, b0): """ inp: t, b0 :float otpt: out, err :float """ stop = False zero = self._mantissa[0] * 0 #one = zero + 1 #1 of type of zero if t < 0: return zero, zero max_a = zero ex = self._exponent if ex >= 0: raise ValueError( 'PN exponent = {} > 0 - inverese Laplce transform does not exist.'.format(ex)) ce1 = -ex-1 pi_ = pi if (t == 0): if (b0 != 0): return zero, zero # b0==0 and t==0: if (ex+2 >= 0): out = self._mantissa[ex+2] err = zero else: out = zero err = zero else: # t > 0 y = [zero,zero] # ring-array of 2 last values[0..1] y[0] = exp(-(b0*b0)/(4*t)) / sqrt(pi_*t) y[1] = erfc( b0 / sqrt(2*t) / sqrt(2+zero) ) for k in range(1, ce1+1): y[(k+1) % 2] = (2*t* y[(k-1) % 2] - b0 * y[k % 2]) / k k = ce1 out = self._mantissa[k-ce1] * y[k % 2] abs_a = [0,0,0] #ring-array of 3 last abs_ak[0..2] abs_a[(ce1+1 -1) % 3] = 0 abs_a[(ce1+1 +0) % 3] = abs(out) dk_maxLocal = -1 while 1: k += 1 y[(k+1) % 2] = (2*t* y[(k-1) % 2] - b0 * y[k % 2]) / k ak = self._mantissa[k-ce1] * y[k % 2] abs_ak = abs(ak) if (abs_ak > max_a): max_a = abs_ak abs_a[(k+1) % 3] = abs_ak if (abs_a[(k+3-1) % 3] < abs_a[k % 3]) \ and (abs_a[k % 3] >= abs_a[(k+1) % 3]) : # local max is here if dk_maxLocal < 0: #start maxLocal = abs_a[k % 3] k_maxLocal = k dk_maxLocal = 0 else: q_maxLocal = abs_a[k % 3]/maxLocal maxLocal = abs_a[k % 3] dk_maxLocal = k-k_maxLocal k_maxLocal = k stop = ( (dk_maxLocal > 0) and (q_maxLocal < 1) and (maxLocal < abs(out)*epsilonPNdig) ) \ or \ ( (k > k_maxLocal + self._max_N // 4) and (abs_ak < abs(out)*epsilonPNdig) ) #last line - for not oscillatory convergent case out = out + ak if (k >= self._max_N-1 + ce1) or stop: break if ( (k <= k_maxLocal+self._max_N-1 // 4) and (dk_maxLocal > 0) and (q_maxLocal<1) ): err = maxLocal * q_maxLocal **((k-k_maxLocal)/dk_maxLocal) else: err = abs_ak err = err + (max_a+out) * epsilonPNdig return out, err #invTr05exp_b0_LaplPN() #### def const_Qu(self, n, a): """ Qu_n(a) ========== returns PN constants connected with Bessel functions I and K; k-th digit (k=0,1,...) = Γ(n+k+(1/2)) / (Γ(n-k+(1/2)) * k! * (2*a) **k) n, a - float, although in practice n=0 and n=1 -- not yet tested """ _4nn = 4*n*n mant = self._mantissa[:] mant[0] = mant[0]*0 + 1 # 1.0 of type of mant[0] for k in range(1, self._max_N): mant[k] = mant[k-1] * (_4nn -(2*k-1)*(2*k-1)) / (8*a*k) return PolyNum(mant, 0) #$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ def getMantissaExponent_fromStr(s, sep, max_N): ''' max_N is used to generate constatnts '~,~' is equivalent to '~,' Some test are also near __init__ ''' if not s: return [], 0 me = s.replace(" ", "").split('**') if len(me) > 2: raise ValueError( "{!r} error - only 1 exponent allowed here.".format(s)) if len(me) > 1: # expo = int(me[1].strip('()')) expo = me[1] expo = expo.lstrip('(').rstrip(')') expo = int(expo) # -2 else: expo = 0 mLR_ = me[0].split('*')[0] # '(~1.2~,2.~-0.3~)' - ignore '*(~1~0~)' # const if mLR_[0] == 'c': mant = [] if mLR_ == 'const:(~2~,-4~4~-4~4~...~)': for __ in range(max_N//2): mant.extend([4,-4]) mant[0] = 2 elif mLR_ == 'const:(~1~,2~2~2~2~...~)': mant = [1] + [2 for __ in range(max_N-1)] else: raise ValueError( "{!r} error - unknow PN const format.".format(s)) return mant[:max_N], expo # values mLR_ = mLR_.replace(sep+',','\`,') # '(~1.2\`,2.~-0.3~)' (for case '(\`,2.~-0.3~)') # mLR = mLR_.strip('()'+sep).split('\`,') # ['1.2\`, '2.~-0.3'] mLR = mLR_.lstrip('('+sep).rstrip(')'+sep) mLR = mLR.split('\`,') # ['1.2\`, '2.~-0.3'] if not mLR[0]: #ex. for input '\`,2.~-0.3' -> mLR[0] == '' mL=['0'] else: mL = mLR[0].split(sep) expo += len(mL) - 1 #ex. for input '(~-1.~2.2~,-3.3~)' -> expo += 1 if len(mLR) > 1: # mLR_.strip(sep) '~,~' case -> '~... mR = mLR[1].lstrip(sep) mR = mR.split(sep) # ['2.', '-0.3'] else: mR = [] m = mL + mR # ['1.2', '2.', '-0.3'] # mant = [float(d) for d in m] mant = [flt(d) for d in m] return mant, expo ### def erfc(x): if x == 0.0: return 1.0 xPosit = ( x>0 ) x = abs(x) max_I = 2 / epsilonPNdig #big float sqrt_pi = sqrt(pi) if x < 1.5: wk0 = exp(-x*x) * 2 * x / sqrt_pi out1 = wk0 k = 1 while 1: wk0 = wk0 * 2 * x*x / (2*k+1) stop = (abs(wk0) * max_I < abs(out1)) out1 = out1 + wk0 k += 1 if stop: #or (k+2 > max_N) break erfc0 = abs(1-out1) #err0 = abs(wk0) else: wk0 = exp(-x*x) / (sqrt_pi*x) out = wk0 abs_wk0 = abs(wk0) k = 1 while 1: wk0 = wk0 * (2*k-1)/(-2*x*x) stop = (abs(wk0) > abs_wk0) if stop: #neutralize last added item: out = out - wk0 * (-2*x*x)/(2*k-1) else: abs_wk0 = abs(wk0) stop = (abs_wk0 * max_I < abs(out)) out = out + wk0 k += 1 if stop: #or (k+2 > max_N) break erfc0 = abs(out) #err111 = abs_wk0 return erfc0 if xPosit else 2 - erfc0 #################### MantPN functions ################################# #################### CONST ############################################ _powersOfTwo = [1] while _powersOfTwo[-1] < PolyNumConf_max_N: _powersOfTwo += [2 * _powersOfTwo[-1]] _powersOfTwo = tuple(_powersOfTwo) # (1, 2, 4, 8, 16, 32, 64, 128) ####################################################################### def mantPN_val(mant, x): """ Evaluate (Horner's scheme) a MantPN ''mant'' at specific value ''x'' - scalar or other objects (like PN). Type of x is most important. """ ### >>> p(PolyNum([5.])) # 3. * 5.**0 + 0 * 5.**(-1)1 + 1 * 5.**(-2) + 2 * 5.**(-3) ### PolyNum('(~3.056~)') # # not important case, consuming long time x_1 = 1/x y = mant[-1] # mantR = mant[::-1] #transcrypt doesn''t like ''[::-1]'' mantR = list(reversed(mant)) for p1 in mantR[1:]: y = y * x_1 + p1 return y def mantPN_mul(x, h, N): """ Find the product of two mantPN. Parameters ---------- x, h : mantPN objects, both must have N elements! Returns ------- out : (x * h)[:N] """ zero = (x[0]+h[0])*0 # zero of common-type of x and y y = [zero for __ in range(N)] # y = [zero]*N for k in range(N): for j in range(k+1): y[k] = y[k] + x[k-j]*h[j] return y #y = np.convolve(x, h, mode='full') #'same' - max(M, N); 'full' - M+N-1 #return y[:N] def mantPN_inv(x, N): """ Find an inversion of mantPN. Parameters ---------- x : mantPN object, must have N elements! Returns ------- y : (1 / x)[:N] """ zero = x[0] * 0 LoN = 0 while _powersOfTwo[LoN] < N: LoN += 1 # 2**LoN >= N #print(f'''N={N}, 2** {_powersOfTwo[LoN]} -- {_powersOfTwo}''') v = [zero for __ in range(_powersOfTwo[LoN])] # v = [zero] * _powersOfTwo[LoN] y = v[:] if len(x) < len(v): # case if N != _powersOfTwo[LoN] x = (x + v)[:len(v)] y[0] = 1 / x[0] #x[0]**(-1), despite of x[0] type for w in range(LoN): nw = _powersOfTwo[w] for k in range(nw): v[k] = zero # unnecessary? for j in range(nw): v[k] = v[k] + y[j] * x[k+nw-j] for k in range(nw): y[k+nw] = zero # unnecessary? for j in range(k+1): y[k+nw] = y[k+nw] - y[k-j] * v[j] #? if abs(y[k+nw]) >= sqrt(MaxFloat): y[k+nw] = sqrt(MaxFloat) ? return y[:N] def mantPN_sqrt(x, N): """ sqrt(x)[:N] """ if x[0] <= 0: raise ValueError("{sqrt(x)} 1st digit of x is not positive: {}".format(x[0])) zero = x[0] * 0 LoN = 1 while _powersOfTwo[LoN] < N: LoN += 1 # 2**LoN >= N y = [zero for __ in range(_powersOfTwo[LoN])] # [zero] * _powersOfTwo[LoN] odwr = y[:] outp = y[:] if len(x) < len(y): # case if N != _powersOfTwo[LoN] x = (x + y)[:len(y)] outp[0] = sqrt(x[0]) for w in range(1, LoN+1): odwr = mantPN_inv(outp, _powersOfTwo[w]) y = mantPN_mul(odwr, x, _powersOfTwo[w]) #mantPN_div(outp,x,_powersOfTwo[w]) for j in range(_powersOfTwo[w-1]): outp[j] = (outp[j] + y[j]) / 2 # 1/2 despite of () type for j in range(_powersOfTwo[w-1], _powersOfTwo[w]): outp[j] = y[j] / 2 return outp[:N] def mantPN_power_real(x, a, N): """ self **a ======== where exists self[0] **a (a is rather float, fract or int) """ zero = x[0] * 0 outp = [zero for __ in range(N)] # y = [zero]*N outp[0] = x[0] **a a += 1 for k in range(1, N): for j in range(1, k+1): outp[k] = outp[k] + x[j] * outp[k-j] * (a*j/k - 1) outp[k] = outp[k] / x[0] return outp def mantPN_exp(x, N): """ """ zero = x[0]*0 exp_x0 = exp(x[0]) outp = [exp_x0 * x1 for x1 in x] outp[0] = exp_x0 for k in range(1, N): #outp[k] = exp_x0 * x[k] for j in range(1, k): outp[k] += x[k-j] * outp[j] * (1 - (zero + j)/k) return outp def mantPN_ln(x, N): """ ln(x) ===== """ zero = x[0]*0 ln_x0 = log(x[0]) outp = x[:] outp[0] = ln_x0 for k in range(1, N): #outp[k] = x[k] for j in range(1, k): outp[k] -= x[j] * outp[k-j] * (1 - (zero + j)/k) outp[k] = outp[k] / x[0] return outp #$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ \`\`\``