; BIGNUM 0.18 - October 3, 2004 ; by Michael Kuyumcu in NetLogo 2.1beta2 ; Begin: September 24, 2004 ; *********************************************** ; CALCULATION ROUTINES ; *********************************************** to-report imag_part [a] ; IN : a large complex number (or number string) ; OUT: its real part ifelse not list? a [report "0"] ; if a is a real number, trivial [report last a] end to-report real_part [a] ; IN : a large complex number (or number string) ; OUT: its real part ifelse not list? a [if number? a [set a number_to_string a] ; trivial case report a] ; if a is a real number, also trivial [report first a] end to-report real_nth_power [a nth] ; IN : a large real number (or number string) a and a parameter ; nth for the (integer) power a is to be raised to ; OUT: a raised to the nth power locals [result] ; Norming the numbers and handling trivial cases if string? nth [set nth read-from-string nth] ; This must be changed if nth is allowed ; to be a bignum someday (as soon as new ; logarithms are available if not integer? nth [report "real_nth_power: n must be an integer, but is not."] if negative? nth [report "real_nth_power: n must not be negative, but is."] if zero? nth [report "1"] if one? nth [report a] if number? a [set a number_to_string a] ; This must be changed if nth is allowed ; to be a bignum someday (as soon as new ; logarithms are available set a sans_zeroes a if zero? a [report "0"] if one? a [report "1"] if two? a [report real_mul a a] ;output-print "real_nth_power: a: " + a + ", nth: " + nth report real_nth_power_produce a nth end to-report real_nth_power_produce [a nth] ; IN : a large real number (or number string) a and a parameter ; nth for the (integer) power a is to be raised to ; ! : this function calls itself recursively locals [counter times result rest timcnt] if zero? nth [report "1"] if one? nth [report a] set times int ((ln nth) / (ln 2)) ; This must be changed if nth is allowed ; to be a bignum someday (as soon as new ; logarithms are available ;output-print "real_nth_power_produce: times: " + times set counter "0" set result a set timcnt "1" while [less? counter times] [set result real_mul result result set counter int_add counter "1" set timcnt int_add timcnt timcnt ] set rest int_sub nth timcnt set rest read-from-string rest ; This must be changed if nth is allowed ; to be a bignum someday (as soon as new ; logarithms are available ;output-print "real_nth_power_produce: rest: " + rest set result real_mul result real_nth_power_produce a rest report result end to-report real_cube_root [a dec_places] ; IN : a large real number (or number string) a and a parameter ; specifying the desired accuracy of the result (decimal places) ; OUT: the cube (3rd) root of a with the desired accuracy report real_nth_root a 3 dec_places ; the cube root is just a special case end to-report real_sqrt [a dec_places] ; IN : a large real number (or number string) a and a parameter ; specifying the desired accuracy of the result (decimal places) ; OUT: the square (2nd) root of a with the desired accuracy report real_nth_root a 2 dec_places ; the square root is just a special case end to-report real_nth_root [a nth dec_places] ; IN : one string (or number) containing a large decimal number and one string (or number) ; parameter nth defining the maximum number of wanted decimal places of the result ; OUT: one string containing the nth (n an integer) root of as calculated with ; a Newton-Rhapson iteration (qudratic-convergence approximation) ; A first approximation is arrived at by dividing a by 2. ; ! reports "" (the empty string) by an illegal division is attempted (by 0) locals [result nth_sub1 approx_old approx_new done? h1 h2 h3 h4 minus_nth] ; The norming of the numbers if string? dec_places [set dec_places read-from-string dec_places] if number? a [set a number_to_string a] set a sans_zeroes a if negative? a [report "real_nth_root: cannot draw the real root of a negative number"] if zero? a [report "0"] if one? a [report "1"] if number? nth [set nth number_to_string nth] set nth sans_zeroes nth if not integer? nth [report "real_nth_root: n must be an integer."] ifelse negative? nth [set minus_nth true set nth butfirst nth] [set minus_nth false] ; Special (trivial) cases if zero? nth [report "1"] if one? nth [ifelse minus_nth [report real_div "1" a dec_places] [report a]] set nth_sub1 int_sub nth "1" ; Now for the initial approximation (guess) set approx_old real_div a nth (dec_places + 2) ; Now we have an approximation for the Newton-Rhapson method that isn't too far away ; from a value needed for the procedure to converge properly (approx_old) set approx_new "0" set done? false while [not done?] ; while the required accuracy has not been reached... [set h1 real_mul nth_sub1 approx_old set h2 real_nth_power approx_old nth_sub1 set h3 real_div a h2 (dec_places + 4) set h4 real_add h1 h3 set approx_new real_div h4 nth (dec_places + 4) set done? ((equal? approx_old approx_new) or (decimal_sameness approx_old approx_new dec_places >= dec_places)) set approx_old approx_new] ; In approx_new we now have an approximation of 1/b of the desired accuracy ifelse minus_nth [set result real_div "1" approx_new dec_places] [set result approx_new] set result sans_zeroes (rounded result dec_places) report result end to-report complex_modulus_squared [a] ; IN : one 2-element list of a complex numbers a ; OUT: the (real-valued) squared modulus of a locals [a_real a_imag result] ifelse not list? a [if number? a [set a number_to_string a] ; trivial case report a] ; if a is a string, also trivial [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] set result real_add (real_mul a_real a_real) (real_mul a_imag a_imag) report result end to-report complex_modulus [a dec_places] ; IN : one 2-element list of a complex numbers a ; OUT: the (real-valued) modulus of a locals [a_real a_imag square_sum result] if string? dec_places [set dec_places read-from-string dec_places] ifelse not list? a [if number? a [set a number_to_string a] ; trivial case report a] ; if a is a string, also trivial [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] set square_sum real_add (real_mul a_real a_real) (real_mul a_imag a_imag) set result real_sqrt square_sum dec_places report result end to-report complex_conjugate [a] ; IN : one 2-element list of a complex numbers a ; OUT: the complex conjugate of a locals [a_real a_imag] ifelse not list? a [if number? a [set a number_to_string a] ; a seems to be a real number set a_real a set a_imag "0"] ; now a is a string [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] report list a_real (change_sign a_imag) end to-report complex_div [a b dec_places] ; IN : two lists of real and imaginary parts of complex numbers in strings ; or two strings of real numbers ; OUT: the complex quotient of both numbers as a 2-element list of strings locals [a_real a_imag b_real b_imag result_real result_imag det] ifelse not list? a [if number? a [set a number_to_string a] ; a seems to be a real number set a_real a set a_imag "0"] [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] ifelse not list? b [if number? b [set b number_to_string b] ; b seems to be a real number set b_real b set b_imag "0"] [set b_real first b set b_imag last b ; b is a complex number if number? b_real [set b_real number_to_string b_real] if number? b_imag [set b_imag number_to_string b_imag]] ; Special trivial cases if (zero? b_real) and (zero? b_imag) [report "complex_div: Division by complex zero not implemented yet :)"] if (one? b_real) and (zero? b_imag) [report list a_real a_imag] ;output-print "complex_div before rounding: a_real: " + a_real +", a_imag: " + a_imag + ", b_real: " + b_real + ", b_imag: " + b_imag ; Normal cases. First, round arguments set a_real rounded a_real dec_places set a_imag rounded a_imag dec_places set b_real rounded b_real dec_places set b_imag rounded b_imag dec_places ;output-print "complex_div: a_real: " + a_real +", a_imag: " + a_imag + ", b_real: " + b_real + ", b_imag: " + b_imag set det real_add (real_mul b_real b_real) (real_mul b_imag b_imag) set result_real real_add (real_mul a_real b_real) (real_mul a_imag b_imag) set result_real real_div result_real det dec_places set result_imag real_sub (real_mul b_real a_imag) (real_mul a_real b_imag) set result_imag real_div result_imag det dec_places report list result_real result_imag end to-report complex_mul [a b] ; IN : two lists of real and imaginary parts of complex numbers in strings ; or two strings of real numbers ; OUT: the complex product of both numbers as a 2-element list of strings locals [a_real a_imag b_real b_imag result_real result_imag] ifelse not list? a [if number? a [set a number_to_string a] ; a seems to be a real number set a_real a set a_imag "0"] [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] ifelse not list? b [if number? b [set b number_to_string b] ; b seems to be a real number set b_real b set b_imag "0"] [set b_real first b set b_imag last b ; b is a complex number if number? b_real [set b_real number_to_string b_real] if number? b_imag [set b_imag number_to_string b_imag]] set result_real real_sub (real_mul a_real b_real) (real_mul a_imag b_imag) set result_imag real_add (real_mul a_imag b_real) (real_mul a_real b_imag) report list result_real result_imag end to-report complex_add [a b] ; IN : two lists of real and imaginary parts of complex numbers in strings ; or two strings of real numbers ; OUT: the complex sum of both numbers as a 2-element list of strings locals [a_real a_imag b_real b_imag result_real result_imag] ifelse not list? a [if number? a [set a number_to_string a] ; a seems to be a real number set a_real a set a_imag "0"] [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] ifelse not list? b [if number? b [set b number_to_string b] ; b seems to be a real number set b_real b set b_imag "0"] [set b_real first b set b_imag last b ; b is a complex number if number? b_real [set b_real number_to_string b_real] if number? b_imag [set b_imag number_to_string b_imag]] set result_real real_add a_real b_real set result_imag real_add a_imag b_imag report list result_real result_imag end to-report complex_sub [a b] ; IN : two lists of real and imaginary parts of complex numbers in strings ; or two strings of real numbers ; OUT: the complex difference of both numbers as a 2-element list of strings locals [a_real a_imag b_real b_imag result_real result_imag] ifelse not list? a [if number? a [set a number_to_string a] ; a seems to be a real number set a_real a set a_imag "0"] [set a_real first a set a_imag last a ; a is a complex number if number? a_real [set a_real number_to_string a_real] if number? a_imag [set a_imag number_to_string a_imag]] ifelse not list? b [if number? b [set b number_to_string b] ; b seems to be a real number set b_real b set b_imag "0"] [set b_real first b set b_imag last b ; b is a complex number if number? b_real [set b_real number_to_string b_real] if number? b_imag [set b_imag number_to_string b_imag]] set result_real real_sub a_real b_real set result_imag real_sub a_imag b_imag report list result_real result_imag end ; ------------------------------- above written October 3, 2004 ------------------------------- to-report real_add [a b] ; IN : two big signed or unsigned real numbers or number strings a and b ; OUT: their sum (a + b) as a string ; HOW: the adding is performed as a subtraction of the negative version of b ; This way, only one routine has to be written (real_sub) locals [minus_a minus_b swap int_a int_b counter frc_a frc_b pos summ dig_a dig_b len_a len_b sign max_len result result_frc result_int carry] ;output-print " real_add: a: "+ a + ", b: " + b if number? a [set a number_to_string a] if number? b [set b number_to_string b] ; Special trivial cases if a = change_sign b [report "0"] if zero? a [report b] if zero? b [report a] if (integer? a) and (integer? b) [report int_add (int_part a) (int_part b)] ; Sign matters set minus_a negative? a set minus_b negative? b if minus_a [set a butfirst a] if minus_b [set b butfirst b] ; Kill leading zeroes in a and b set a int_sans_leading_zeroes a set b int_sans_leading_zeroes b if minus_a and not minus_b [report real_sub b a] ; if only one of the numbers is negative, return the subtraction result if minus_b and not minus_a [report real_sub a b] ; Calculation of the sum to commence. set int_a int_part a set int_b int_part b set frc_a fract_part a if frc_a = "0" [set frc_a "0.0"] set frc_a butfirst frc_a set frc_a butfirst frc_a ; omit leading "0." set frc_b fract_part b if frc_b = "0" [set frc_b "0.0"] set frc_b butfirst frc_b set frc_b butfirst frc_b ; omit leading "0." set len_a length frc_a set len_b length frc_b set max_len max list len_a len_b ;output-print "real_add: a: " + a + ", b: " + b + ", minus_a: " + minus_a + ", minus_b: " + minus_b ;output-print " int_a: " + int_a + ", frc_a: " + frc_a + ", int_b: " + int_b + ", frc_b: " + frc_b set result_frc "" set carry 0 set counter 0 while [counter < max_len] [set pos max_len - counter ; Which digits are to be subtracted? ifelse pos > len_a [set dig_a 0 ] [set dig_a read-from-string substring frc_a (pos - 1) pos] ifelse pos > len_b [set dig_b 0 ] [set dig_b read-from-string substring frc_b (pos - 1) pos] set summ dig_a + dig_b + carry ifelse summ > 9 [set summ summ - 10 set carry 1] [set carry 0 ] set result_frc (summ + result_frc) set counter counter + 1] set result_int int_add int_a carry set result_int int_add result_int int_b set result result_int + "." + result_frc ifelse (minus_a and minus_b) [set result "-" + result] [set result result] ; both a and b were positive ;if substring result 0 2 = "--" [set result butfirst result set result butfirst result] report result end to-report real_sub [a b] ; IN : two big signed or unsigned real numbers or number strings a and b ; OUT: their difference (a - b) as a string ; HOW: the adding is performed as a subtraction of the negative version of b ; This way, only one routine has to be written (real_sub) locals [minus_a minus_b swap int_a int_b counter frc_a frc_b pos diff dig_a dig_b len_a len_b sign max_len result result_frc result_int carry] if number? a [set a number_to_string a] if number? b [set b number_to_string b] ; Special trivial cases if equal? a b [report "0"] if zero? a [report change_sign b] if zero? b [report a] if (integer? a) and (integer? b) [report int_sub (int_part a) (int_part b)] ; Sign matters set minus_a (negative? a) set minus_b (negative? b) if minus_a [set a butfirst a] if minus_b [set b butfirst b] ; Kill leading zeroes in a and b set a int_sans_leading_zeroes a set b int_sans_leading_zeroes b if (not minus_a) and minus_b [report real_add a b] ; this is basically an addition if minus_a and (not minus_b) [report "-" + real_add a b] ; addition deeper into the negative realm if minus_a and minus_b [set swap a set a b set b swap] ; subtract b from a ; Standardize the order of the arguments! set len_a length a set len_b length b set sign "" ; If a < b, swap the numbers and remember in SIGN that the result will have to be negated if greater? b a [set swap a set a b set b swap set sign "-"] ; Calculation of difference. We know that the absolute value of b is equal or less than that of a, ; so b is usually the shorter string (or, at most, of equal length) set int_a int_part a set int_b int_part b set frc_a fract_part a set frc_b fract_part b if frc_a = "0" [set frc_a "0.0"] if frc_b = "0" [set frc_b "0.0"] set frc_a butfirst frc_a set frc_a butfirst frc_a ; omit leading "0." set frc_b butfirst frc_b set frc_b butfirst frc_b ; omit leading "0." set len_a length frc_a set len_b length frc_b set max_len max list len_a len_b ;output-print "real_sub: a: " + a + ", b: " + b + ", minus_a: " + minus_a + ", minus_b: " + minus_b ;output-print " int_a: " + int_a + ", frc_a: " + frc_a + ", int_b: " + int_b + ", frc_b: " + frc_b set result_frc "" set carry 0 set counter 0 while [counter < max_len] [set pos max_len - counter ; Which digits are to be subtracted? ifelse pos > len_a [set dig_a 0 ] [set dig_a read-from-string substring frc_a (pos - 1) pos] ifelse pos > len_b [set dig_b 0 ] [set dig_b read-from-string substring frc_b (pos - 1) pos] set diff dig_a - dig_b - carry ifelse diff < 0 [set diff diff + 10 set carry 1] [set carry 0 ] set result_frc (diff + result_frc) set counter counter + 1] set result_int int_sub int_a carry set result_int int_sub result_int int_b set result sign + result_int + "." + result_frc if substring result 0 2 = "--" [set result butfirst result set result butfirst result] report result end to-report real_div [a b dec_places] ; IN : two strings (or numbers) containing large decimal numbers and one parameter ; defining the maximum number of wanted decimal places of the result ; OUT: one string containing the quotient of a and b as calculated with ; the Newton-Rhapson procedure (qudratic-convergence approximation) ; A first approcximation is arrived at by a built-in division of a few ; first significant digits of b. a/b is written as aá1/b, and 1/b is ; calculated first. As a rather close first guess for 1/b is crucial for the ; convergence of Newton's method, we examine b closely. ; ! reports "" (the empty string) by an illegal division is attempted (by 0) locals [result dec_places_b b_int l dec_pos_b dec_pos_a result_sign approx_old approx_new approx_help minus_a minus_b counter done? found] if number? a [set a number_to_string a] if number? b [set b number_to_string b] set a sans_zeroes a set b sans_zeroes b ; Negative values? ifelse negative? a [set minus_a true set a butfirst a] [set minus_a false] ifelse negative? b [set minus_b true set b butfirst b] [set minus_b false] ifelse ((minus_a and not minus_b) or (minus_b and not minus_a)) [set result_sign "-"] [set result_sign "" ] ; Special cases - special treatment if zero? b [report "Division by zero not implemented."] ; Illegal attempt at division by zero if zero? a [report "0"] ; 0/something = 0 if one? b [report result_sign + a] ; Division by one (trivial case) if (b = "0.5" or b = ".5") [report result_sign + real_mul a 2] ; Division by one half if (b = "0.25" or b = ".25") [report result_sign + real_mul a 4] ; Division by one quarter if (b = "0.1" or b = ".1") [report result_sign + power10 a 1] ; Division by one tenth ; Remove decimal point from denominator set dec_pos_b position "." b ifelse dec_pos_b != false [set b_int int_part b] [set b_int b] ; b was an integer already ifelse b_int != "0" [set approx_help "" + (10 / read-from-string substring b_int 0 1) set l (length b_int) - 1 if approx_help = "10" [set l l - 1 set approx_help "1"] set dec_pos_b position "." approx_help if dec_pos_b != false [set approx_help remove-item dec_pos_b approx_help] ; omit decimal point set counter 0 while [counter < l] [set approx_help "0" + approx_help set counter counter + 1] set approx_old "0." + approx_help ;output-print "real_div: 1st approx.: " + approx_help + ", b_int: " + b_int + ", l: " + l + ", approx_old: " + approx_old ] ; Now we have a fractional value for b below 1 (0.xyz...) ; and must find the first digit after the decimal point <> 0 [set found false set counter 0 while [(counter < length b) and not found] [ifelse (substring b (dec_pos_b + counter + 1) (dec_pos_b + counter + 2)) = "0" [set counter counter + 1] [set found true]] set approx_help substring b (dec_pos_b + counter + 1) (dec_pos_b + counter + 2) set approx_help "" + (10 / read-from-string approx_help) ifelse counter > 0 [set approx_old real_tp10 approx_help counter] [set approx_old approx_help] ;output-print "real_div: b_int zero: counter: " + counter + ", approx_help: " + approx_help + ", approx_old: " + approx_old ] ; Now we have an approximation for the Newton-Rhapson method that isn't too far away ; from a value needed for the procedure to converge properly (approx_old) set approx_new "0" set done? false while [not done?] ; while the required accuracy has not been reached... [set approx_new rounded (real_mul approx_old b) (dec_places + 2) set approx_new real_sub "2" approx_new set approx_new rounded (real_mul approx_new approx_old) (dec_places + 2) set done? ((equal? approx_old approx_new) or (decimal_sameness approx_old approx_new dec_places >= dec_places)) set approx_old approx_new] ; In approx_new we now have an approximation of 1/b of the desired accuracy set result result_sign + (real_mul a approx_new) set result sans_zeroes (rounded result dec_places) report result end ; ------------------------------- above written October 2, 2004 ------------------------------- to-report rounded [a places] ; IN: one strings (or number) containing large decimal number and a ; second number (#places) specifying to how many decimal digits a ; should be rounded. ; OUT: a string with the rounded value of a locals [a_frc dig_a sign_a] if number? a [set a number_to_string a] if empty? a [report "rounded: cannot round an empty argument"] ifelse negative? a [set sign_a "-" set a butfirst a] [set sign_a ""] set a_frc fract_part a if a_frc = "0" [set a_frc "0.0"] set a_frc butfirst a_frc set a_frc butfirst a_frc if length a_frc < (places + 1) [report sign_a + a] set a_frc substring a_frc 0 (places + 1) set dig_a last a_frc ifelse dig_a < "5" [report sign_a + (int_part a) + "." + (butlast a_frc)] [set a (int_part a) + "." + (butlast a_frc) ;output-print "rounded: a: " + a report sign_a + (real_add a (power10 1 (- places)))] end to-report decimal_sameness [a b places] ; IN: two strings (or numbers) containing large decimal numbers and a ; third number (#places) specifying how many decimal digits a and b ; should have in common. ; OUT: an integer telling the number of corresponding digits in a row ; In the case of places being 0, only the integer part is considered ; In this case, the output gives TRUE (the same) or FALSE (not the same) locals [a_int a_frc b_int b_frc same_places char len] if number? a [set a number_to_string a] if number? b [set b number_to_string b] if string? places [set places read-from-string places] set a sans_zeroes a set b sans_zeroes b set a_int int_part a set b_int int_part b set a_frc fract_part a if a_frc = "0" [set a_frc "0.0"] set b_frc fract_part b if b_frc = "0" [set b_frc "0.0"] ;output-print "decimal_sameness: a_int: " + a_int + ", b_int: " + b_int + ", a_frc: " + a_frc + ", b_frc: " + b_frc ifelse places = 0 [report a_int = b_int] ; only integer parts are relevant [if a_int != b_int [report 0] set len min list length a_frc length b_frc set same_places 0 set char 2 while [(same_places < places) and (char < len)] [if (substring a_frc char (char + 1)) = (substring b_frc char (char + 1)) [set same_places same_places + 1] set char char + 1 ] report same_places ] end to-report real_mul [a b] ; IN: two strings (or numbers) containing large decimal numbers ; OUT: one string containing the PRODUCT of the numbers ; HOW: The number strings are treated as integers, and after the calculation ; of the integer product, the decimal point is inserted locals [dec_pos_a dec_pos_b dec_places_a dec_places_b dec_places_allt result res_dec_pos minus_a minus_b] if number? a [set a number_to_string a] if number? b [set b number_to_string b] if (zero? a) or (zero? b) [report "0"] ; Negative values? Make them positive and remember signs ifelse negative? a [set minus_a true set a butfirst a] [set minus_a false] set a sans_zeroes a ifelse negative? b [set minus_b true set b butfirst b] [set minus_b false] set b sans_zeroes b ;output-print "a new: " + a + ", minus_a: " + minus_a + ", b new: " + b + ", minus_b: " + minus_b ; Is there a decimal point? set dec_pos_a position "." a ; Position of the decimal points (if any) set dec_pos_b position "." b ; If there are decimal places, remember how many and delete the decimal point ; That is, we pass integer values to the main multiplication routine and ; re-insert the decimal point afterwards. ifelse dec_pos_a = false [set dec_places_a 0] [set a real_sans_trailing_zeroes a set dec_pos_a position "." a set dec_places_a length a - dec_pos_a - 1 set a remove-item dec_pos_a a] ifelse dec_pos_b = false [set dec_places_b 0] [set b real_sans_trailing_zeroes b set dec_pos_b position "." b set dec_places_b length b - dec_pos_b - 1 set b remove-item dec_pos_b b] ; How many decimal places are there altogether? set dec_places_allt dec_places_a + dec_places_b set result int_mul a b ; multiply the numbers as if they were integer numbers while [length result < dec_places_allt] [set result "0" + result] if dec_places_allt > 0 [set res_dec_pos (length result) - dec_places_allt set result (substring result 0 res_dec_pos) + "." + substring result res_dec_pos (length result)] set result real_sans_trailing_zeroes result if zero? result [report "0"] ifelse (minus_a and not minus_b) or (minus_b and not minus_a) [report "-" + result] [report result] end to-report power10 [number n] ; IN : a number or number string, a positive or negative power of 10 (n) ; OUT: this number multiplied with the respective power of 10 locals [len_number len_diff dec_pos_old dec_pos_new zeroes sign] if number? number [set number number_to_string number] ; if numerical, make it into a string if string? n [set n read-from-string n] ; if not numerical, make it into a number set number sans_zeroes number if n = 0 [report number] ; case n = 0: no change to the number handed ifelse negative? number [set sign "-" set number butfirst number] [set sign ""] ifelse n > 0 [report sign + real_tp10 number n] ; n positive or zero [set len_number length int_part number set len_diff len_number + n - 1 ; "+n" is, in fact, a subtraction of n, as n is negative if len_diff <= 0 [set zeroes string "0" (abs len_diff) set number zeroes + number] ; prepend a number of zeroes for moving of the decimal point set dec_pos_old position "." number ifelse dec_pos_old = false [set dec_pos_new length number + n] ; case of integer number [set number remove-item dec_pos_old number ; decimal case set dec_pos_new dec_pos_old + n] set number (substring number 0 dec_pos_new) + "." + (substring number dec_pos_new length number) set number sans_zeroes number report sign + number ; decimal number ] end to-report string [str places] ; IN : a string to be duplicated #places times. Places is a number ; OUT: the duplicated string locals [result counter] if number? str [set str number_to_string str] if string? places [set places read-from-string places] if places < 0 [report "Cannot produce negative multiples of a string."] if places = 0 [report ""] if places = 1 [report str] report string_produce str places end to-report string_produce [s places] ; IN : a string to be duplicated #places times. Places is a number ; OUT: the multiplied string ; ! : recursive doubling of s locals [counter result times rest] if number? s [set s number_to_string s] if places = 1 [report s] set result s set times int (ln places / ln 2) ; how many duplications do we need? set counter 0 while [counter < times] [set result result + result set counter counter + 1] set rest places - length result if rest > 0 [set result (result + string_produce s rest)] ; recursive call for the rest report result end to-report real_tp10 [number n] ; MULTIPLY A REAL OR INTEGER NUMBER WITH A POWER OF 10 ; IN: one string (or number) containing a large (signed/unsigned) REAL number ; OUT: one string containing the product of the number and the n-th power of ten (n >= 0) locals [dec_pos result] if number? number [set number number_to_string number] ; if numerical, make it into a string if string? n [set n read-from-string n] ; if not numerical, make it into a number set dec_pos position "." number ifelse dec_pos = false [set result int_tp10 number n] ; integer special case (no decimal point manipulation) [set result remove-item dec_pos number set result int_tp10 result n set dec_pos dec_pos + n set result (substring result 0 dec_pos) + "." + (substring result dec_pos length result)] set result real_sans_leading_zeroes result set result real_sans_trailing_zeroes result report result end to-report int_tp10 [number n] ; MULTIPLY A NUMBER WITH A POWER OF 10 ; IN: one string (or number) containing a large (signed/unsigned) INTEGER number, one number n (power) ; OUT: one string containing the product of the number and the n-th power of ten (n >= 0) locals [add_on dec_pos counter len_num rest_dig] if number? number [set number number_to_string number] ; if numerical, make it into a string if string? n [set n read-from-string n] ; if not numerical, make it into a number if n < 0 [report "int_tp10 cannot work with a negative number of zeroes"] report number + string "0" n end to-report int_sub [a b] ; IN: two strings or numbers containing large (positive or negative) integer numbers a, b ; OUT: one string containing a large integer number c, the DIFFERENCE of a and b (= a - b) locals [minus_a minus_b result swap sign diff carry digit dig_a dig_b pos_a pos_b len_a len_b maxlen counter] if is-number? a [set a number_to_string a] ; this routine can be called with arguments of type string or number if is-number? b [set b number_to_string b] ; if numbers are handed, they are converted into strings first if zero? a [report b] if zero? b [report a] set minus_a negative? a set minus_b negative? b if minus_a [set a butfirst a] if minus_b [set b butfirst b] ; Kill leading zeroes in a and b set a int_sans_leading_zeroes a set b int_sans_leading_zeroes b if equal? a b [report "0"] if (not minus_a) and minus_b [report int_add a b] ; simple addition if minus_a and (not minus_b) [report "-" + int_add a b] ; addition into the negative realm if minus_a and minus_b [set swap a set a b set b swap] ; subtract b from a set len_a length a set len_b length b set sign "" ; If a < b, swap the numbers and remember in SIGN that the result will have to be negated if less? a b [set swap a set a b set b swap set sign "-"] ; Calculation of difference. We know that the absolute value of b is equal or less than that of a, ; so b is usually the shorter string (or, at most, of equal length) set result "" set carry 0 set len_a length a ; neccessary to compute again since the order of a and b set len_b length b ; might have been swapped set counter 0 while [counter < len_a] [set pos_a len_a - counter - 1 set pos_b len_b - counter - 1 set dig_a read-from-string (substring a pos_a (pos_a + 1)) ifelse pos_b >= 0 [set dig_b read-from-string (substring b pos_b (pos_b + 1))] [set dig_b 0] set diff dig_a - dig_b - carry ifelse diff < 0 [set carry 1 set diff 10 + diff] [set carry 0] if (pos_a != 0) or (diff != 0) [set result (diff + result)] set counter counter + 1] report sign + result end to-report int_add [a b] ; IN: two strings or numbers containing large (positive or negative) integer numbers a, b ; OUT: one string containing a large integer number c, the SUM of a and b locals [minus_a minus_b result summ carry digit dig_a dig_b pos_a pos_b len_a len_b maxlen counter] if is-number? a [set a number_to_string a] ; this routine can be called with string or number arguments if is-number? b [set b number_to_string b] ; if numbers are handed, they are converted into strings first if zero? a [report b] if zero? b [report a] set minus_a negative? a set minus_b negative? b if minus_a [set a butfirst a] if minus_b [set b butfirst b] ; Kill leading zeroes in a and b set a int_sans_leading_zeroes a set b int_sans_leading_zeroes b if minus_a and not minus_b [report int_sub b a] ; if only one of the numbers is negative, return the subtraction result if minus_b and not minus_a [report int_sub a b] ; MAIN CALCULATION LOOP (adding with carry) set result "" set carry 0 set len_a length a set len_b length b set maxlen max list len_a len_b set counter 0 while [counter < maxlen] [set pos_a len_a - counter - 1 set pos_b len_b - counter - 1 ifelse pos_a >= 0 [set dig_a read-from-string (substring a pos_a (pos_a + 1))] [set dig_a 0] ifelse pos_b >= 0 [set dig_b read-from-string (substring b pos_b (pos_b + 1))] [set dig_b 0] set summ dig_a + dig_b + carry ifelse summ > 9 [set carry 1 set summ summ - 10] [set carry 0] set result (summ + result) set counter counter + 1] if carry = 1 [set result ("1" + result)] ;output-print "Result int_add: " + result ifelse (minus_a and minus_b) [report "-" + result] [report result] ; both a and b were positive end to-report int_mul [a b] ; IN: two strings OR numbers containing large (positive or negative) integer numbers a, b ; OUT: one STRING containing a large integer number c, the PRODUCT of a and b locals [minus_a minus_b sign ] ; the signs of the numbers. The numbers will be multiplied without the signs first if is-number? a [set a number_to_string a] ; this routine can be called with string or number arguments if is-number? b [set b number_to_string b] ; if numbers are handed, they are converted into strings first set minus_a ((substring a 0 1) = "-") set minus_b ((substring b 0 1) = "-") if minus_a [set a butfirst a] if minus_b [set b butfirst b] set a sans_zeroes a set b sans_zeroes b if zero? a [report "0"] if one? a [report b] if zero? b [report "0"] if one? b [report a] ;output-print "int_mul: a: " + a + ", b: "+ b ifelse ((minus_a and minus_b) or (not minus_a and not minus_b)) [set sign ""] ; positive result (no sign) [set sign "-"] ; negative result ifelse (length a < 4) and (length b < 4) [report sign + (read-from-string a) * (read-from-string b)] ; simple case, small numbers [report sign + int_mult_toom_cook a b] end to-report int_mult_toom_cook [a b] ; IN: a and b are strings containing large positive numbers ; OUT: the product of a and b as a string ; HOW: Each integer's digits are divided into 3 equal length numbers, ; w digits wide. For example: a = 37865512 will be divided to ; 037,865,512, w = 3. The numbers can now be represented as ; a = (A1, A2, A3) and b = (B1, B2, B3). The result will be: ; c = aáb = (C1,C2,C3,C4,C5) while ; C1 = A1áB1 ; C2 = A1áB2+A2áB1 ; C3 = A1áB3+A2áB2+A3áB1 ; C4 = A2áB3+A3áB2 ; C5 = A3áB3 locals [result len_a len_b len_parts_a len_parts_b a1 a2 a3 b1 b2 b3 c1 c2 c3 c4 c5 a10? a20? a30? b10? b20? b30? ; Flags: Is a1 zero? Is a2 zero? etc. a11? a21? a31? b11? b21? b31? ; Flags: Is a1 one? Is a2 one? etc. sum_1 sum_2 sum_3 counter] if zero? a [report "0"] if zero? b [report "0"] if one? a [report b] if one? b [report a] if two? a [report int_add b b] if two? b [report int_add a a] if number? a [set a number_to_string a] if number? b [set b number_to_string b] ; Divide up a into three chunks of same size ; But first, make the numbers same size. set len_a length a set len_b length b if (len_b > len_a) [set len_parts_a len_b - len_a set counter 0 while [counter < len_parts_a] [set a "0" + a set counter (counter + 1)] ] if (len_a > len_b) [set len_parts_b len_a - len_b set counter 0 while [counter < len_parts_b] [set b "0" + b set counter (counter + 1)] ] ; Now we can split a and b up into three chunks of the same size set len_a length a ifelse len_a > 3 [set len_parts_a int (len_a / 3) if (len_a mod 3) > 0 [set len_parts_a 1 + len_parts_a] ;output-print "LŠnge von a: " + len_a + ", LŠnge der Teile: " + len_parts set a3 substring a (len_a - len_parts_a) len_a set a2 substring a (len_a - 2 * len_parts_a) (len_a - len_parts_a) set a1 substring a 0 min list len_parts_a (len_a - 2 * len_parts_a) while [length a1 < len_parts_a] [set a1 "0" + a1]] [set a3 a set a2 "0" set a1 "0"] ; Divide up b into three chunks of same size, too set len_b length b ifelse len_b > 3 [set len_parts_b int (len_b / 3) if (len_b mod 3) > 0 [set len_parts_b 1 + len_parts_b] ;output-print "LŠnge von b: " + len_b + ", LŠnge der Teile: " + len_parts set b3 substring b (len_b - len_parts_b) len_b set b2 substring b (len_b - 2 * len_parts_b) (len_b - len_parts_b) set b1 substring b 0 min list len_parts_b (len_b - 2 * len_parts_b) while [length b1 < len_parts_b] [set b1 "0" + b1]] [set b3 b set b2 "0" set b1 "0"] ; FIRST CASE: ; Here we have small numbers that can be multiplied with the internal * operator ifelse (len_parts_a <= 1) and (len_parts_b <= 1) ; numbers are small enough to be handled by [set a1 read-from-string a1 ; internal multiplication operator * set a2 read-from-string a2 ; So peel them out of their string shell set a3 read-from-string a3 ; and make them ready for numerical set b1 read-from-string b1 ; multiplication set b2 read-from-string b2 set b3 read-from-string b3 set c1 (a1 * b1) ; These c1 to c5 values are the constituents set c2 (a1 * b2 + a2 * b1) ; of the final product which we will embed set c3 (a1 * b3 + a2 * b2 + a3 * b1) ; again into a string jacket. set c4 (a2 * b3 + a3 * b2) set c5 (a3 * b3) set result int_mul c1 int_tp10 1 (4 * len_parts_a) set result result + int_mul c2 int_tp10 1 (3 * len_parts_a) set result result + int_mul c3 int_tp10 1 (2 * len_parts_a) set result result + int_mul c4 int_tp10 1 len_parts_a set result result + c5 ;output-print "c1: " + c1 + ", c2: " + c2 + ", c3: " + c3 + ", c4: " + c4 + ", c5: " + c5 report "" + result] ; SECOND CASE: ; Now the numbers are too big to be handled by the built-in multiplication operator *, ; so the routine makes a recursive call to itself with smaller factors [set a10? zero? a1 set b10? zero? b1 set a11? one? a1 set b11? one? b1 set a20? zero? a2 set b20? zero? b2 set a21? one? a2 set b21? one? b2 set a30? zero? a3 set b30? zero? b3 set a31? one? a3 set b31? one? b3 ; Calculate C1 and take into cosideration special cases like one factor is zero or one ifelse (a10? or b10?) [set c1 "0"] [ifelse a11? [set c1 b1] [ifelse b11? [set c1 a1] [set c1 int_mult_toom_cook a1 b1] ] ] ; Calculate C5 in the same manner ifelse (a30? or b30?) [set c5 "0"] [ifelse a31? [set c5 b3] [ifelse b31? [set c5 a3] [set c5 int_mult_toom_cook a3 b3] ] ] ; Calculate C2 in a similar manner, but try to simplify even more cases ; Split c2 up into two summed parts, sum_1 and sum_2. First consider sum_1 ifelse (a10? or b20?) [set sum_1 "0"] [ifelse a11? [set sum_1 b2] [ifelse b21? [set sum_1 a1] [set sum_1 int_mult_toom_cook a1 b2] ] ] ifelse (a20? or b10?) [set sum_2 "0"] [ifelse a21? [set sum_2 b1] [ifelse b11? [set sum_2 a2] [set sum_2 int_mult_toom_cook a2 b1] ] ] ifelse (zero? sum_1) [set c2 sum_2] [ifelse zero? sum_2 [set c2 sum_1] [set c2 int_add sum_1 sum_2] ] ; Calculate C4 in the same manner, and again try to simplify even more cases ; Split c4 up into two summed parts, sum_1 and sum_2. First consider sum_1 ifelse (a20? or b30?) [set sum_1 "0"] [ifelse a21? [set sum_1 b3] [ifelse b31? [set sum_1 a2] [set sum_1 int_mult_toom_cook a2 b3] ] ] ifelse (a30? or b20?) [set sum_2 "0"] [ifelse a31? [set sum_2 b2] [ifelse b21? [set sum_2 a3] [set sum_2 int_mult_toom_cook a3 b2] ] ] ifelse (zero? sum_1) [set c4 sum_2] [ifelse zero? sum_2 [set c4 sum_1] [set c4 int_add sum_1 sum_2] ] ; Finally, calculate the more complex C3 case, and again try to simplify as many trivial cases as possible ; Split c3 up into three summed parts, sum_1, sum_2, and sum_3. First consider sum_1 ifelse (a10? or b30?) [set sum_1 "0"] [ifelse a11? [set sum_1 b3] [ifelse b31? [set sum_1 a1] [set sum_1 int_mult_toom_cook a1 b3] ] ] ifelse (a20? or b20?) [set sum_2 "0"] [ifelse a21? [set sum_2 b2] [ifelse b21? [set sum_2 a2] [set sum_2 int_mult_toom_cook a2 b2] ] ] ifelse (a30? or b10?) [set sum_3 "0"] [ifelse a31? [set sum_3 b1] [ifelse b11? [set sum_3 a3] [set sum_3 int_mult_toom_cook a3 b1] ] ] set c3 int_add (int_add sum_1 sum_2) sum_3 ; Now calculate the final result from the partial sums set c1 int_tp10 c1 (4 * len_parts_a) ; multiplies c1 and 10 to the power of (4 * len_parts_a) set c2 int_tp10 c2 (3 * len_parts_a) set c3 int_tp10 c3 (2 * len_parts_a) set c4 int_tp10 c4 len_parts_a set result int_add c1 c2 set result int_add result c3 set result int_add result c4 set result int_add result c5 report result ] end to-report negative? [n] if number? n [set n number_to_string n] if empty? n [report false] report (substring n 0 1) = "-" end to-report positive? [n] if number? n [set n number_to_string n] if empty? n [report false] set n sans_zeroes n ifelse substring n 0 1 = "-" [report false] [set n sans_zeroes n report not zero? n] end to-report integer? [n] if number? n [set n number_to_string n] if empty? n [report false] set n sans_zeroes n report ((position "." n) = false) end to-report int_part [n] ; IN : a large signed or unsigned number or number string n ; OUT: the INTeger part of that number (no rounding) locals [dec_pos_n sign_n] ; Format norming if number? n [set n number_to_string n] if empty? n [report ""] ifelse negative? n [set sign_n "-" set n butfirst n] [set sign_n ""] set n real_sans_leading_zeroes n if substring n 0 1 = "." [set n "0" + n] set dec_pos_n position "." n ifelse dec_pos_n = false [report sign_n + n] ; number was an integer already [report sign_n + (substring n 0 dec_pos_n)] end to-report fract_part [n] ; IN : a large signed or unsigned number or number string n ; OUT: the fractional part of that number (no rounding) locals [dec_pos_n result] if number? n [set n number_to_string n] if empty? n [report ""] set dec_pos_n position "." n ifelse dec_pos_n = false [report "0"] ; number was an integer already [set result substring n dec_pos_n length n ifelse zero? result [report "0"] [report "0" + result] ] end ; *********************************************** ; CALCULATION HELPER ROUTINES ; *********************************************** to-report number_to_string [n] ; IN : a NetLogo number ; OUT: this number's representation as a string ; ! : NetLogo codes some numbers in "scientific notation", ; for instance "1.255242E-6" locals [e_pos pow num] set n "" + n ; make it into a string set e_pos position "E" n ifelse e_pos = False [report n] ; not in scientific notation [set pow substring n (e_pos + 1) (length n) set num substring n 0 e_pos report power10 num pow ] end to-report zero? [n] ; this routine returns FALSE (n is not zero) or TRUE (it is zero) locals [c l] if number? n [report (n = 0)] if empty? n [report false] if substring n 0 1 = "-" [set n butfirst n] ; discar minus sign ; Decimal number? Ignore the decimal point set l position "." n if l != false [set n remove-item l n] set c 0 set l length n while [c < l] [if (substring n c (c + 1)) != "0" [report false] set c c + 1] report true end to-report one? [n] ; this routine returns FALSE (n is not 1) or TRUE (it is 1) if number? n [report (n = 1)] if empty? n [report false] set n sans_zeroes n report n = "1" end to-report two? [n] ; this routine returns FALSE (n is not 2) or TRUE (it is 2) locals [c l] if number? n [report (n = 2)] if empty? n [report false] set n sans_zeroes n report n = "2" end to-report int_sans_leading_zeroes [n] ; Eliminate leading zeroes in an integer number (string) if number? n [report number_to_string n] if empty? n [report ""] while [(substring n 0 1 = "0") AND (length n > 1)] [set n butfirst n] report n end to-report real_sans_trailing_zeroes [n] locals [pos] ; removes the trailing zeroes from a number of type real ; it is important that the number handed to this routine really does contain a decimal point! if number? n [report number_to_string n] if empty? n [report ""] set pos position "." n if pos != false [while [(substring n (length n - 1) (length n)) = "0"] [set n but-last n] if position "." n = (length n - 1) [set n butlast n] ] report n end to-report real_sans_leading_zeroes [n] locals [dec_pos sign result] if number? n [set n number_to_string n] if empty? n [report ""] ifelse (position "-" n) = 0 [set sign "-" set n butfirst n] [set sign ""] ; Is there a decimal point in the number? set dec_pos position "." n if (dec_pos = 0) [report sign + "0" + n] ifelse dec_pos != false [if dec_pos > 0 [set result (int_sans_leading_zeroes substring n 0 dec_pos) + "." + substring n (dec_pos + 1) length n]] [set result int_sans_leading_zeroes n] if (zero? result) or (result = "") [report "0"] report sign + result end to-report sans_zeroes [n] ; Takes any number (string or integer or real) and omits leading and trailing zeroes. ; IN : number or number string (any format) ; OUT: that number with trailing and leading zeroes removed if number? n [set n number_to_string n] if empty? n [report ""] set n real_sans_leading_zeroes n set n real_sans_trailing_zeroes n ifelse zero? n [report "0"] [report n] end ; NUMBER VALUE COMPARISON OPERATORS to-report greater? [a b] ; IN : two (integer or real) numbers or number strings a, b ; OUT: a boolean telling whether (TRUE) or not (FALSE) a is strictly greater than b locals [a_int b_int int_greater fract_greater a_frc b_frc swap len_a len_b sign_a sign_b] if number? a [set a number_to_string a] if number? b [set b number_to_string b] set a sans_zeroes a set b sans_zeroes b ; Negative values? ifelse negative? a [set sign_a "-" set a butfirst a] [set sign_a ""] ifelse negative? b [set sign_b "-" set b butfirst b] [set sign_b ""] ; Simple cases; if both numbers are negative, swap them if sign_a = "" and sign_b = "-" [report true] if sign_a = "-" and sign_b = "" [report false] if sign_a = "-" and sign_b = "-" [set swap a set a b set b swap set sign_a "" set sign_b ""] ; this line only in case it is needed later set a_int int_part a set len_a length a_int set a_frc fract_part a ; in case the fractional pieces decide set b_int int_part b set len_b length b_int set b_frc fract_part b set int_greater (len_a > len_b) or ((len_a = len_b) and (a_int > b_int)) set fract_greater (a_frc > b_frc) ;output-print "greater_than: a: " + a + ", b: " + b + ", a_int: " + a_int + ", b_int: " + b_int + ", len_a: " + len_a + ", len_b: " + len_b ;output-print " a_frc: " + a_frc + ", b_frc: " + b_frc + ", int_greater: " + int_greater + ", fract_greater: " + fract_greater report int_greater or ((equal? a_int b_int) and fract_greater) end to-report less? [a b] report greater? b a end to-report less_or_equal? [a b] report (greater? b a) or (equal? a b) end to-report greater_or_equal? [a b] ; IN : two numbers or number strings a, b ; OUT: a boolean telling whether (TRUE) or not (FALSE) a is strictly greater than b report (greater? a b) or (equal? a b) end to-report equal? [a b] ; IN : two numbers or number strings a, b ; OUT: a boolean telling whether (TRUE) or not (FALSE) a is equal to b if number? a [set a number_to_string a] if number? b [set b number_to_string b] set a sans_zeroes a set b sans_zeroes b report (a = b) end to-report real_abs [a] ; IN : one number or number string containing a large signed or unsigned real or integer number ; OUT: the number with positve (i. e. no) sign if number? a [set a number_to_string a] set a sans_zeroes a ifelse negative? a [report butfirst a] [report a] end to-report real_sgn [a] ; IN : one number or number string containing a large signed or unsigned real or integer number ; OUT: the sign of the number ("1" positive, "0" for zero, and "-1" for a negative number) if number? a [set a number_to_string a] set a sans_zeroes a if negative? a [report "-1"] ifelse zero? a [report "0"] [report "1"] end to-report change_sign [a] ; IN : one number or number string containing a large signed or unsigned real or integer number ; OUT: the same number with its sign reversed ( + -> -, - -> +) if number? a [set a number_to_string a] if empty? a [report "change_sign cannot deal with an empty argument."] if zero? a [report "0"] ifelse positive? a [report "-" + a] [report butfirst a] end ; *********************************************************** ; * TEST SUITE ; *********************************************************** to test_suite ; Test locals [t1 millisec counter a b c] clear-output output-print "BIGNUM Test Suite" ; Long float Addition set a rounded "101.92882827271827272662626262515144295673535272" precisio set b rounded "5.982772726262827267255242423636635352452422423095947345325" precisio output-print "\n1000 long float additions (" + precisio + " digits precision): " set t1 timer set counter 0 while [counter < 1000] [set c real_add a b set counter counter + 1] set millisec (timer - t1) output-print "Result of " + a + " + " + b + "=" output-print (rounded c precisio) + "\n1000 calculations needed " + millisec + " seconds." ; Long float Subtraction output-print "\n1000 long float subtractions (" + precisio + " digits precision): " set a rounded "101.92882827271827272662626262515144295673535272" precisio set b rounded "5.982772726262827267255242423636635352452422423095947345325" precisio set t1 timer set counter 0 while [counter < 1000] [set c real_sub a b set counter counter + 1] set millisec (timer - t1) output-print "Result of " + a + " - " + b + "=" output-print (rounded c precisio) + "\n1000 calculations needed " + millisec + " seconds." ; Long float Multiplication output-print "\n10 long float multiplications (" + precisio + " digits precision): " set a rounded "101.92882827271827272662626262515144295673535272" precisio set b rounded "5.982772726262827267255242423636635352452422423095947345325" precisio set t1 timer set counter 0 while [counter < 10] [set c real_mul a b set counter counter + 1 ] set millisec (timer - t1) output-print "Result of " + a + " * " + b + "=" output-print (rounded c precisio) + "\n10 calculations needed " + millisec + " seconds." ; Long float Division output-print "\n10 long float divisions (" + precisio + " digits precision): " set a rounded "101.92882827271827272662626262515144295673535272" precisio set b rounded "5.982772726262827267255242423636635352452422423095947345325" precisio set t1 timer set counter 0 while [counter < 10] [set c real_div a b precisio set counter counter + 1 output-type counter + " "] set millisec (timer - t1) output-print "\nResult of " + a + " / " + b + "=" output-print (rounded c precisio) + "\n10 calculations needed " + millisec + " seconds." ; Long integer Addition set b "10192882827271827272662626262515144295673535272" set a "5982772726262827267255242423636635352452422423095947345325" output-print "\n1000 long integer additions (full precision): " set t1 timer set counter 0 while [counter < 1000] [set c int_add a b set counter counter + 1] set millisec (timer - t1) output-print "Result of " + a + " + " + b + "=" output-print (rounded c precisio) + "\n1000 calculations needed " + millisec + " seconds." ; Long integer Subtraction set b "10192882827271827272662626262515144295673535272" set a "5982772726262827267255242423636635352452422423095947345325" output-print "\n1000 long integer subtractions (full pecision): " set t1 timer set counter 0 while [counter < 1000] [set c int_sub a b set counter counter + 1] set millisec (timer - t1) output-print "Result of " + a + " - " + b + "=" output-print (rounded c precisio) + "\n1000 calculations needed " + millisec + " seconds." ; Long integer Multplication set b "1019288282727182729568877773535272" set a "598277272626282355144145325" output-print "\n10 long integer multiplications (full precision): " set t1 timer set counter 0 while [counter < 10] [set c int_mul a b set counter counter + 1] set millisec (timer - t1) output-print "Result of " + a + " * " + b + "=" output-print (rounded c precisio) + "\n10 calculations needed " + millisec + " seconds." end @#$#@#$#@ GRAPHICS-WINDOW 1010 453 1185 649 4 4 18.333333333333332 1 10 1 1 1 0 CC-WINDOW 5 671 1203 766 Command Center OUTPUT 17 49 1004 646 TEXTBOX 17 17 167 35 Tests BUTTON 1013 100 1186 133 Test BIGNUM routines test_suite NIL 1 T OBSERVER T NIL SLIDER 1013 51 1185 84 Precisio Precisio 1 201 21 2 1 NIL @#$#@#$#@ WHAT IS IT? ----------- This is a collection of reporters for doing calculations with arbitrary-precision integer and floating-point real and complex numbers. In the collection you find the following reporters which accept both STRINGS and normal NetLogo numbers as arguments. In the case of complex numbers, for each complex number, one 2-element list of strings (real part and imaginary part) are passed. Results are ALWAYS passed back to the caller as STRINGS (real floating-point numbers) or 2-element lists of strings (complex numbers). Example: set c real_mul "1288272.827727272727" "188272727.292929299229922992292971" (Result for c: "242546638773587.705228117835542710057003926002101917") or, one example for complex numbers: set c complex_mul ["127.18872" "-19.277271882"] ["-0.2377161" "8.7272662"] (Result for c: "[138.0030768415969884 1114.5923349676927002]") int_add a b : adds two large - signed or unsigned - integer numbers int_sub a b : subtracts b from a (signed or unsigned) int_mul a b : multiplies large - signed or unsigned - numbers real_add a b : adds two large - unsigned or signed - floating-point numbers real_sub a b : subtracts a large floating-point number a from b real_mul a b : multiplies two large floating-point numbers real_div a b precision : divides two large floating-point numbers with #precision digits real_nth_power a nth : raises a (a real value) to the non-negative integer power n real_nth_root a nth accuracy : calculates the nth (n an integer) root of a : which is correct to #accuracy decimal digits real_sqrt a accuracy : calculates the square root of a : which is correct to #accuracy decimal digits real_cube_root a accuracy : calculates the cube (3rd) root of real number a complex_add a b: sum of two complex numbers (handed as 2-element lists of strings), for instance: set c complex_add ["12.11" "-0.1"] ["-9.881" "2.112"] The first list element is the real part, the 2nd the imaginary complex_sub a b: calculates the complex difference a - b. complex_mul a b: calculates the comples product a * b complex_div a b precision: calculates a / b (b must not equal ["0" "0"]) to the desired accuracy. complex_conjugate a: computes the complex conjugate of a complex number a complex_modulus_squared a: computes the squared modulus of a complex number a complex_modulus a precision: gives the square root of complex_modulus_squared : with the desired precision real_part a : gives the real part of the complex number a imag_part a : gives the imaginary part of the complex number a real_sgn a : reports the sign of a number (1 positive, 0 zero, -1 negative) change_sign a : multiplies a and (-1) real_abs a : returns -a for negative a, and a for positive a int_part a : reports the integer portion of a floating point number a fract_part a : reports the fractional portion of a rounded a dec_places : reports a rounded to #dec_places decimal places integer? a : reports TRUE if a has no decimal point, FALSE otherwise positive? a : reports TRUE if a > 0, FALSE otherwise negative? a ; reports TRUE if a < 0, FALSE otherwise zero? a : reports TRUE if a equals 0, FALSE otherwise equal? a b : reports TRUE if a equals b, FALSE otherwise greater? a b : reports TRUE if the value of a is greater than b's less? a b : reports TRUE if a < b, FALSE otherwise greater_or_equal? a b : analogous less_or_equal? a b : analogous string text times : reports a string with #text appended #times times. : e.g. string "H" 10 -> "HHHHHHHHHH" power10 number times : reports the product of number and 10 raised to the : the power of #times decimal_sameness a b dec_places : for #dec_places > 0, reports the number of decimal places shared by a and b (continuous) : for #dec_places = 0, reports whether a and b have the same integer value : This is reporter is useful for iterations HOW IT WORKS ------------ Pass arguments to the reporters mentioned above as normal NetLogo numbers or as STRINGS (for arbitrary-length and -precision numbers) containing the numbers to be processed in base 10 notation. The most elaborate routines are int_mul for the multiplication of large integers (which is used for floating-point numbers as well). That reporter uses the TOOM-COOK algorithm which desiplays a much better runtime behavior than multiplying with the longhand "school" method, i.e. O(n^1.485) vs. O(n^2). At the heart of the floating point division algorithm is a standard Newton-Rhapson iteration. As a sufficiently close first approximation of the solution for 1/b in (a/b = a * 1/b) is crucial for this algorithm, the routines tries to find out a close guess, but it appears to be clumsy. Rewrite it with a better approach - and mail it to me! :) info@noemanetz.de Thanks! My motivation in writing these routines was my interest in fractal geometry. For producing my own fractal exploration software, I need customized routines that can handle numbers with arbitrary precision (length). Most of the programs that have been written so far for creating Mandelbrot sets, Julia sets and images of Newton-iterated roots of complex functions are not capable of handling really high-precision numbers. As a result, the user can zoom into the calculated fractal image only a few times before the built-in arithmetic becomes too "coarse" for deep detail. But realizing how slow NetLogo is when it comes to string manipulation, I have decided to only use it as a rapid prototyping environment and will rewrite the routines in a different programming language after testing of the routines. HOW TO USE IT ------------- Run the test suite by pressing the button on the interface tab. You may use the routines provided here as you like. THINGS TO NOTICE ---------------- Normally, when doing iterations and checking whether the old approximated value is sufficiently close to the new one, one subtracts both, calculated |Xold-Xnew| and checks if this is less a certain threshold epsilon. With the reporter decimal_sameness, the check is made "graphically" by comparing the significant decimal places instead of a full-blown subtraction and comparison (which in itself is another comparison). If you are interested in this approach, have a look at the WHILE loop in the real_div reporter. THINGS TO TRY ------------- call real_mul with really large numbers. Time the result. See how slow NetLogo string handling is. EXTENDING THE MODEL ------------------- Write reporters that draw roots (square, cubic, arbitrary). Implement reporters for arbitrary powers (fractional), logarithms, trigonometric functions and a few for number theory like Li(n). Implement representations for continued fractions. Add routines for complex number calculations. NETLOGO FEATURES ---------------- - RELATED MODELS -------------- - CREDITS AND REFERENCES ---------------------- - @#$#@#$#@ default true 0 Polygon -7566196 true true 150 5 40 250 150 205 260 250 ant true 0 Polygon -7566196 true true 136 61 129 46 144 30 119 45 124 60 114 82 97 37 132 10 93 36 111 84 127 105 172 105 189 84 208 35 171 11 202 35 204 37 186 82 177 60 180 44 159 32 170 44 165 60 Polygon -7566196 true true 150 95 135 103 139 117 125 149 137 180 135 196 150 204 166 195 161 180 174 150 158 116 164 102 Polygon -7566196 true true 149 186 128 197 114 232 134 270 149 282 166 270 185 232 171 195 149 186 Polygon -7566196 true true 225 66 230 107 159 122 161 127 234 111 236 106 Polygon -7566196 true true 78 58 99 116 139 123 137 128 95 119 Polygon -7566196 true true 48 103 90 147 129 147 130 151 86 151 Polygon -7566196 true true 65 224 92 171 134 160 135 164 95 175 Polygon -7566196 true true 235 222 210 170 163 162 161 166 208 174 Polygon -7566196 true true 249 107 211 147 168 147 168 150 213 150 arrow true 0 Polygon -7566196 true true 150 0 0 150 105 150 105 293 195 293 195 150 300 150 bee true 0 Polygon -256 true false 152 149 77 163 67 195 67 211 74 234 85 252 100 264 116 276 134 286 151 300 167 285 182 278 206 260 220 242 226 218 226 195 222 166 Polygon -16777216 true false 150 149 128 151 114 151 98 145 80 122 80 103 81 83 95 67 117 58 141 54 151 53 177 55 195 66 207 82 211 94 211 116 204 139 189 149 171 152 Polygon -7566196 true true 151 54 119 59 96 60 81 50 78 39 87 25 103 18 115 23 121 13 150 1 180 14 189 23 197 17 210 19 222 30 222 44 212 57 192 58 Polygon -16777216 true false 70 185 74 171 223 172 224 186 Polygon -16777216 true false 67 211 71 226 224 226 225 211 67 211 Polygon -16777216 true false 91 257 106 269 195 269 211 255 Line -1 false 144 100 70 87 Line -1 false 70 87 45 87 Line -1 false 45 86 26 97 Line -1 false 26 96 22 115 Line -1 false 22 115 25 130 Line -1 false 26 131 37 141 Line -1 false 37 141 55 144 Line -1 false 55 143 143 101 Line -1 false 141 100 227 138 Line -1 false 227 138 241 137 Line -1 false 241 137 249 129 Line -1 false 249 129 254 110 Line -1 false 253 108 248 97 Line -1 false 249 95 235 82 Line -1 false 235 82 144 100 bird1 false 0 Polygon -7566196 true true 2 6 2 39 270 298 297 298 299 271 187 160 279 75 276 22 100 67 31 0 bird2 false 0 Polygon -7566196 true true 2 4 33 4 298 270 298 298 272 298 155 184 117 289 61 295 61 105 0 43 boat1 false 0 Polygon -1 true false 63 162 90 207 223 207 290 162 Rectangle -6524078 true false 150 32 157 162 Polygon -16776961 true false 150 34 131 49 145 47 147 48 149 49 Polygon -7566196 true true 158 33 230 157 182 150 169 151 157 156 Polygon -7566196 true true 149 55 88 143 103 139 111 136 117 139 126 145 130 147 139 147 146 146 149 55 boat2 false 0 Polygon -1 true false 63 162 90 207 223 207 290 162 Rectangle -6524078 true false 150 32 157 162 Polygon -16776961 true false 150 34 131 49 145 47 147 48 149 49 Polygon -7566196 true true 157 54 175 79 174 96 185 102 178 112 194 124 196 131 190 139 192 146 211 151 216 154 157 154 Polygon -7566196 true true 150 74 146 91 139 99 143 114 141 123 137 126 131 129 132 139 142 136 126 142 119 147 148 147 boat3 false 0 Polygon -1 true false 63 162 90 207 223 207 290 162 Rectangle -6524078 true false 150 32 157 162 Polygon -16776961 true false 150 34 131 49 145 47 147 48 149 49 Polygon -7566196 true true 158 37 172 45 188 59 202 79 217 109 220 130 218 147 204 156 158 156 161 142 170 123 170 102 169 88 165 62 Polygon -7566196 true true 149 66 142 78 139 96 141 111 146 139 148 147 110 147 113 131 118 106 126 71 box true 0 Polygon -7566196 true true 45 255 255 255 255 45 45 45 butterfly1 true 0 Polygon -16777216 true false 151 76 138 91 138 284 150 296 162 286 162 91 Polygon -7566196 true true 164 106 184 79 205 61 236 48 259 53 279 86 287 119 289 158 278 177 256 182 164 181 Polygon -7566196 true true 136 110 119 82 110 71 85 61 59 48 36 56 17 88 6 115 2 147 15 178 134 178 Polygon -7566196 true true 46 181 28 227 50 255 77 273 112 283 135 274 135 180 Polygon -7566196 true true 165 185 254 184 272 224 255 251 236 267 191 283 164 276 Line -7566196 true 167 47 159 82 Line -7566196 true 136 47 145 81 Circle -7566196 true true 165 45 8 Circle -7566196 true true 134 45 6 Circle -7566196 true true 133 44 7 Circle -7566196 true true 133 43 8 circle false 0 Circle -7566196 true true 35 35 230 line true 0 Line -7566196 true 150 0 150 300 person false 0 Circle -7566196 true true 155 20 63 Rectangle -7566196 true true 158 79 217 164 Polygon -7566196 true true 158 81 110 129 131 143 158 109 165 110 Polygon -7566196 true true 216 83 267 123 248 143 215 107 Polygon -7566196 true true 167 163 145 234 183 234 183 163 Polygon -7566196 true true 195 163 195 233 227 233 206 159 sheep false 15 Rectangle -1 true true 90 75 270 225 Circle -1 true true 15 75 150 Rectangle -16777216 true false 81 225 134 286 Rectangle -16777216 true false 180 225 238 285 Circle -16777216 true false 1 88 92 spacecraft true 0 Polygon -7566196 true true 150 0 180 135 255 255 225 240 150 180 75 240 45 255 120 135 thin-arrow true 0 Polygon -7566196 true true 150 0 0 150 120 150 120 293 180 293 180 150 300 150 truck-down false 0 Polygon -7566196 true true 225 30 225 270 120 270 105 210 60 180 45 30 105 60 105 30 Polygon -8716033 true false 195 75 195 120 240 120 240 75 Polygon -8716033 true false 195 225 195 180 240 180 240 225 truck-left false 0 Polygon -7566196 true true 120 135 225 135 225 210 75 210 75 165 105 165 Polygon -8716033 true false 90 210 105 225 120 210 Polygon -8716033 true false 180 210 195 225 210 210 truck-right false 0 Polygon -7566196 true true 180 135 75 135 75 210 225 210 225 165 195 165 Polygon -8716033 true false 210 210 195 225 180 210 Polygon -8716033 true false 120 210 105 225 90 210 turtle true 0 Polygon -7566196 true true 138 75 162 75 165 105 225 105 225 142 195 135 195 187 225 195 225 225 195 217 195 202 105 202 105 217 75 225 75 195 105 187 105 135 75 142 75 105 135 105 wolf false 0 Rectangle -7566196 true true 15 105 105 165 Rectangle -7566196 true true 45 90 105 105 Polygon -7566196 true true 60 90 83 44 104 90 Polygon -16777216 true false 67 90 82 59 97 89 Rectangle -1 true false 48 93 59 105 Rectangle -16777216 true false 51 96 55 101 Rectangle -16777216 true false 0 121 15 135 Rectangle -16777216 true false 15 136 60 151 Polygon -1 true false 15 136 23 149 31 136 Polygon -1 true false 30 151 37 136 43 151 Rectangle -7566196 true true 105 120 263 195 Rectangle -7566196 true true 108 195 259 201 Rectangle -7566196 true true 114 201 252 210 Rectangle -7566196 true true 120 210 243 214 Rectangle -7566196 true true 115 114 255 120 Rectangle -7566196 true true 128 108 248 114 Rectangle -7566196 true true 150 105 225 108 Rectangle -7566196 true true 132 214 155 270 Rectangle -7566196 true true 110 260 132 270 Rectangle -7566196 true true 210 214 232 270 Rectangle -7566196 true true 189 260 210 270 Line -7566196 true 263 127 281 155 Line -7566196 true 281 155 281 192 wolf-left false 3 Polygon -6524078 true true 117 97 91 74 66 74 60 85 36 85 38 92 44 97 62 97 81 117 84 134 92 147 109 152 136 144 174 144 174 103 143 103 134 97 Polygon -6524078 true true 87 80 79 55 76 79 Polygon -6524078 true true 81 75 70 58 73 82 Polygon -6524078 true true 99 131 76 152 76 163 96 182 104 182 109 173 102 167 99 173 87 159 104 140 Polygon -6524078 true true 107 138 107 186 98 190 99 196 112 196 115 190 Polygon -6524078 true true 116 140 114 189 105 137 Rectangle -6524078 true true 109 150 114 192 Rectangle -6524078 true true 111 143 116 191 Polygon -6524078 true true 168 106 184 98 205 98 218 115 218 137 186 164 196 176 195 194 178 195 178 183 188 183 169 164 173 144 Polygon -6524078 true true 207 140 200 163 206 175 207 192 193 189 192 177 198 176 185 150 Polygon -6524078 true true 214 134 203 168 192 148 Polygon -6524078 true true 204 151 203 176 193 148 Polygon -6524078 true true 207 103 221 98 236 101 243 115 243 128 256 142 239 143 233 133 225 115 214 114 wolf-right false 3 Polygon -6524078 true true 170 127 200 93 231 93 237 103 262 103 261 113 253 119 231 119 215 143 213 160 208 173 189 187 169 190 154 190 126 180 106 171 72 171 73 126 122 126 144 123 159 123 Polygon -6524078 true true 201 99 214 69 215 99 Polygon -6524078 true true 207 98 223 71 220 101 Polygon -6524078 true true 184 172 189 234 203 238 203 246 187 247 180 239 171 180 Polygon -6524078 true true 197 174 204 220 218 224 219 234 201 232 195 225 179 179 Polygon -6524078 true true 78 167 95 187 95 208 79 220 92 234 98 235 100 249 81 246 76 241 61 212 65 195 52 170 45 150 44 128 55 121 69 121 81 135 Polygon -6524078 true true 48 143 58 141 Polygon -6524078 true true 46 136 68 137 Polygon -6524078 true true 45 129 35 142 37 159 53 192 47 210 62 238 80 237 Line -16777216 false 74 237 59 213 Line -16777216 false 59 213 59 212 Line -16777216 false 58 211 67 192 Polygon -6524078 true true 38 138 66 149 Polygon -6524078 true true 46 128 33 120 21 118 11 123 3 138 5 160 13 178 9 192 0 199 20 196 25 179 24 161 25 148 45 140 Polygon -6524078 true true 67 122 96 126 63 144 @#$#@#$#@ NetLogo 2.1beta2 @#$#@#$#@ @#$#@#$#@ @#$#@#$#@