| Class | Apfp |
| In: |
lib/rapfp/Apfp.rb
|
| Parent: | Object |
###
File: Apfp.rb
Subject: Class for arbitrary precision floating point.
Author: Dennis J. Darland
Date: April 3, 2007
Copyright (C) 2007 Dennis J. Darland#
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
@mant is mantissa except as a whole number - @exp is exponent (of 10)(negative to produce fractions) - @errmant is mantissa of esimated error (always positive) - @errexp is exponent (of 10) of estimated error - word.
# File lib/rapfp/Apfp.rb, line 43
43: def initialize(mant,expt,errmant,errexpt)
44:
45: @mant = mant
46: @expt = expt
47: @errmant = errmant
48: @errexpt = errexpt
49:
50: end
# File lib/rapfp/Apfp.rb, line 233
233: def *(other)
234: a = self.clone
235: b = other.clone
236: cexpt = a.expt + b.expt
237:
238: errmant1 = a.errmant * b.mant.abs
239: errmant2 = b.errmant * a.mant.abs
240: errexpt1 = a.errexpt + b.expt
241: errexpt2 = b.errexpt + a.expt
242: errshift = errexpt1 - errexpt2
243:
244: if errshift > 0 then
245: errshift.times do errmant1 *= 10 end
246: errexpt1 -= errshift
247: elsif errshift < 0
248: errshift.times do errmant2 *= 10 end
249: errexpt2 -= errshift
250: end
251:
252: Apfp.new(a.mant*b.mant,cexpt,errmant1 + errmant2,errexpt1).norm
253:
254: end
# File lib/rapfp/Apfp.rb, line 581
581: def **(other)
582: a = self.clone
583: unless other.kind_of?(Apfp)
584: b = @@rconst.one
585: if other > 0 then
586: other.times do
587: b = a * b.clone
588: end
589: return b
590: elsif other < 0
591: (-other).times do
592: b = a * b.clone
593: end
594: b = @@rconst.one/b.clone
595: return b
596: end
597: return b
598: else
599: b = other
600: if b == (ap_int(0)) then return ap_int(1) end
601: if b == (ap_int(1)) then return a end
602: if a == (ap_int(0)) then return ap_int(0) end
603: return (a.log*b).exp
604: end
605: end
# File lib/rapfp/Apfp.rb, line 198
198: def +(other)
199: a = self.clone
200: b = other.clone
201: if a.mant == 0 then return b end
202: if b.mant == 0 then return a end
203: shift = a.expt - b.expt
204: if (shift > 0) then
205: shift.times do
206: a = a.manttimes10!
207: end
208: a = a.exptadd!(-shift)
209: elsif (shift < 0) then
210: (-shift).times do
211: b = b.manttimes10!
212: end
213: b = b.exptadd!(shift)
214: end
215:
216: shift2 = a.errexpt - b.errexpt
217: if (shift2 > 0) then
218: shift2.times do a = a.errmanttimes10! end
219: a = a.errexptadd!(-shift2)
220: elsif (shift2 < 0) then
221: (-shift2).times do b = b.errmanttimes10! end
222: b = b.errexptadd!(shift2)
223: end
224:
225: Apfp.new(a.mant+b.mant,a.expt,a.errmant.abs+b.errmant.abs,a.errexpt).norm
226:
227: end
# File lib/rapfp/Apfp.rb, line 267
267: def /(other)
268: a = self.clone
269: b = other.clone
270: amant = a.mant
271: if amant == 0 then return a end
272: (2*NUM_DIGITS).times do amant *= 10 end
273: aexpt = a.expt
274: aexpt -= 2*NUM_DIGITS
275: da = a.makerr
276: db = b.makerr
277: (2*NUM_DIGITS).times do da = da.manttimes10! end
278: da.exptadd!(-2*NUM_DIGITS)
279: dc1 = Apfp.new(da.mant.abs/b.mant.abs,da.expt - b.expt,5,-NUM_DIGITS)
280: dc2 = a*db
281: (2*NUM_DIGITS).times do dc2 = dc2.manttimes10! end
282: dc2.exptadd!(-2*NUM_DIGITS)
283: dc3 = Apfp.new(dc2.mant.abs/b.mant.abs,dc2.expt - b.expt, 5, -NUM_DIGITS)
284: (2*NUM_DIGITS).times do dc3 = dc3.manttimes10! end
285: dc3.exptadd!(-2*NUM_DIGITS)
286: dc4 = Apfp.new(dc3.mant.abs/b.mant.abs,dc3.expt - b.expt, 5, -NUM_DIGITS)
287: dc = dc1 + dc4
288: Apfp.new(amant/b.mant,aexpt - b.expt,dc.mant,dc.expt).norm
289:
290: end
# File lib/rapfp/Apfp.rb, line 292
292: def <(other)
293:
294: (self - other).mant < 0
295:
296: end
very rough & fast
# File lib/rapfp/Apfp.rb, line 524
524: def about_log10
525: b = self.norm
526: return ap_int(b.expt)
527: end
# File lib/rapfp/Apfp.rb, line 360
360: def abs
361: if self < (ap_int(0))
362: self.neg.norm
363: else
364: self.norm
365: end
366: end
# File lib/rapfp/Apfp.rb, line 438
438: def acos
439: a = self.norm
440: if a.abs > (ap_int(1)) then
441: # puts "acos out of range"
442: a.display_val("a")
443: return a
444: end
445: b =(@@rconst.pi/@@rconst.two) - a.asin
446: # puts "acos(" + a.to_s + ") = " + b.to_s
447: return b
448: end
# File lib/rapfp/Apfp.rb, line 406
406: def asin
407: a = self.norm
408: if a.abs > (ap_int(1)) then
409: # puts "asin out of range"
410: a.display_val("a")
411: return a
412: end
413: if a == @@rconst.one then
414: return @@rconst.pi/@@rconst.two
415: end
416: s = ApfpSeriesConst.new("asin",a,@@rconst)
417: b = s.compute
418: # puts "asin(" + a.to_s + ") = " + b.to_s
419: return b
420: end
# File lib/rapfp/Apfp.rb, line 421
421: def asin_basic
422: a = self.norm
423: if a.abs > (ap_int(1)) then
424: # puts "asin_basic out of range"
425: a.display_val("a")
426: return a
427: end
428: s = ApfpSeriesConst.new("asin_basic",a,@@rconst)
429: s.compute
430: end
# File lib/rapfp/Apfp.rb, line 450
450: def atan
451: a = self.clone
452: if a >=(ap_int(0)) && (a*a) <= (ap_int(1))
453: how = 1
454: elsif a < (ap_int(0)) && (a*a) <= (ap_int(1))
455: how = 2
456: elsif a > (ap_int(0))
457: a = ap_int(1)/a.clone
458: how = 3
459: else
460: a = ap_int(1)/a.clone
461: how = 4
462: end
463: b = a.atan2
464: # puts "atan2(" + a.to_s + ") = " + b.to_s + " how = " + how.to_s
465: # return b
466: case how
467: when 1,2
468: return b
469: when 3
470: return @@rconst.pi/@@rconst.two - b
471: when 4
472: return @@rconst.pi/@@rconst.minus_two - b
473: end
474: end
# File lib/rapfp/Apfp.rb, line 476
476: def atan2
477: a = self.norm
478: if a < ap_int(-1) || a > ap_int(1) then
479: # puts "Illegal value in atan2"
480: return @@rconst.zero
481: end
482: s = ApfpSeriesConst.new("atan",a,@@rconst)
483: s.compute
484: end
# File lib/rapfp/Apfp.rb, line 389
389: def cos
390: a = self.norm
391: s = ApfpSeriesConst.new("cos",a,@@rconst)
392: s.compute
393: end
# File lib/rapfp/Apfp.rb, line 631
631: def cosh
632:
633: y1 = self.exp
634: y2 = (self * @@rconst.minus_one).exp
635: (y1+y2) / @@rconst.two
636: end
# File lib/rapfp/Apfp.rb, line 170
170: def display_val(label)
171: if @mant < 0 then
172: mant = -@mant
173: sign = "-";
174: s1 = @mant.to_s.size - 1
175: else
176: mant = @mant
177: sign = ""
178: s1 = @mant.to_s.size
179: end
180: if @errmant < 0 then
181: errmant = -@errmant
182: errsign = "-";
183: s2 = @errmant.to_s.size - 1
184: else
185: errmant = @errmant
186: errsign = ""
187: s2 = @errmant.to_s.size
188: end
189:
190: # puts "#{label} = #{sign}0.#{mant.to_s}e#{(@expt+s1).to_s}+/-#{errsign}0.#{errmant.to_s}e#{(@errexpt+s2).to_s}"
191: end
# File lib/rapfp/Apfp.rb, line 326
326: def eq(other)
327: !(self < other) && !(self > other)
328: end
# File lib/rapfp/Apfp.rb, line 432
432: def erf
433: a = self.norm
434: s = ApfpSeriesConst.new("erf",a,@@rconst)
435: s.compute
436: end
# File lib/rapfp/Apfp.rb, line 165
165: def errexptadd!(other)
166: @errexpt += other
167: self
168: end
# File lib/rapfp/Apfp.rb, line 395
395: def exp
396: a = self.norm
397: s = ApfpSeriesConst.new("exp",a,@@rconst)
398: s.compute
399: end
# File lib/rapfp/Apfp.rb, line 400
400: def exp_err
401: a = self.norm
402: s = ApfpSeriesConst.new("exp_err",a,@@rconst)
403: s.compute
404: end
# File lib/rapfp/Apfp.rb, line 323
323: def gt(other)
324: (self-self.makerr) < (other+other.makerr)
325: end
log base 10 of positive value
# File lib/rapfp/Apfp.rb, line 530
530: def log10
531:
532: if self <= (ap_int(0)) then
533: # puts "log10 out of range"
534: self.display_val("self")
535: return self.abs
536: end
537: a = self.norm
538: if a.mant == 1 then
539: return ap_int(a.expt)
540: end
541: # puts "a.mant = |" + a.mant.to_s + "| a.expt = |" + a.expt.to_s + "|"
542: return ap_int(a.expt+a.mant.to_s.size) + Apfp.new(a.mant,-a.mant.to_s.size,a.errmant,a.errexpt).log_basic/@@rconst.log_e_10
543: end
log base 10 of positive value for purpose of error
# File lib/rapfp/Apfp.rb, line 546
546: def log10_err
547:
548: if self <= (ap_int(0)) then
549: # puts "log10_err out of range"
550: self.display_val("self")
551: return self.abs
552: end
553: a = self.norm
554: # puts "a.mant = |" + a.mant.to_s + "| a.expt = |" + a.expt.to_s + "|"
555: return ap_int(a.expt+a.mant.to_s.size) + Apfp.new(a.mant,-a.mant.to_s.size,a.errmant,a.errexpt).log_basic_err/Apfp.new(2302585093,-9,5,-9)
556: end
log base e of value betwween 0 and one
# File lib/rapfp/Apfp.rb, line 487
487: def log_basic
488: a = self.norm
489: if a <= ap_int(0) || a > ap_int(1) then
490: # puts "log_basic out of range"
491: a.display_val("a")
492: return a
493: end
494: # puts "log basic of " + a.to_s
495: s = ApfpSeriesConst.new("log",a,@@rconst)
496: return s.compute
497: end
log base e of value beteew 0 and 1 for purposes of error calc
# File lib/rapfp/Apfp.rb, line 500
500: def log_basic_err
501: a = self.norm
502: if a <= (ap_int(0)) then
503: # puts "log_basic_err out of range"
504: a.display_val("a")
505: return a
506: end
507: s = ApfpSeriesConst.new("log_err",a,@@rconst)
508: return s.compute
509: end
log base e for purpose of calculating log_e_10 in initialize
# File lib/rapfp/Apfp.rb, line 512
512: def log_basic_init
513: a = self.norm
514: if a <= (ap_int(0)) then
515: # puts "log_basic_init out of range"
516: a.display_val("a")
517: return a
518: end
519: s = ApfpSeriesConst.new("log_init",a,@@rconst)
520: return s.compute
521: end
# File lib/rapfp/Apfp.rb, line 320
320: def lt(other)
321: (self+self.makerr) < (other-other.makerr)
322: end
# File lib/rapfp/Apfp.rb, line 256
256: def makerr
257: Apfp.new(@errmant,@errexpt,5,-NUM_DIGITS).norm
258: end
# File lib/rapfp/Apfp.rb, line 193
193: def neg
194: mant = -1*@mant
195: Apfp.new(mant,@expt,@errmant.abs,@errexpt)
196: end
# File lib/rapfp/Apfp.rb, line 94
94: def norm
95: if @mant == 0 then return Apfp.new(0,1,5,-NUM_DIGITS) end
96: mant = @mant.to_s
97: errmant = @errmant.to_s
98: expt = @expt
99: errexpt = @errexpt
100: rmant = mant.reverse
101: sz = rmant.size
102: if sz > NUM_DIGITS then
103: trim = sz - NUM_DIGITS
104: mre = Regexp.new("^([0-9]){" + trim.to_s + "}")
105: if rmant =~ mre then
106: rmant = "#{$'}"
107: m_cnt = "#{$&}".size
108: else
109: m_cnt = 0
110: end
111: else
112: m_cnt = 0
113: end
114:
115: lz_re = Regexp.new("^0+")
116: if rmant =~ lz_re then
117: rmant = "#{$'}"
118: lz_cnt = "#{$&}".size
119: else
120: lz_cnt = 0
121: end
122: mant = rmant.reverse.to_i
123: expt += (lz_cnt + m_cnt)
124: # now process error
125: rerrmant = errmant.reverse
126: sz = rerrmant.size
127: if sz > NUM_ERR_DIGITS then
128: trim = sz - NUM_ERR_DIGITS
129: mre = Regexp.new("^([0-9]){" + trim.to_s + "}")
130: if rerrmant =~ mre then
131: rerrmant = "#{$'}"
132: m_cnt = "#{$&}".size
133: else
134: m_cnt = 0
135: end
136: lz_re = Regexp.new("^0+")
137: if rerrmant =~ lz_re then
138: rerrmant = "#{$'}"
139: lz_cnt = "#{$&}".size
140: else
141: lz_cnt = 0
142: end
143: errmant = rerrmant.reverse.to_i
144: errexpt += (lz_cnt + m_cnt)
145: else
146: errmant = errmant.to_i
147: end
148: Apfp.new(mant,expt,errmant,errexpt)
149: end
# File lib/rapfp/Apfp.rb, line 368
368: def period(per,phase1,phase2)
369: a = self.clone
370: c = a/per
371: # if (c < (@@rconst.zero))
372: # then
373: # c += per
374: # end
375:
376: d = c.frac
377:
378: d2 = d * per
379: e = d2 + phase2
380: return e
381: end
# File lib/rapfp/Apfp.rb, line 260
260: def seterr!(limit)
261: @errmant = limit.mant
262: @errexpt = limit.expt
263: self.norm
264: end
# File lib/rapfp/Apfp.rb, line 383
383: def sin
384: a = self.norm
385: s = ApfpSeriesConst.new("sin",a,@@rconst)
386: s.compute
387: end
# File lib/rapfp/Apfp.rb, line 624
624: def sinh
625:
626: y1 = self.exp
627: y2 = (self * @@rconst.minus_one).exp
628: (y1-y2) / @@rconst.two
629: end
# File lib/rapfp/Apfp.rb, line 615
615: def sqrt
616: if self == @@rconst.zero then
617: return @@rconst.zero
618: end
619: self ** @@rconst.half
620: end
# File lib/rapfp/Apfp.rb, line 607
607: def sqrt_err
608: if self == Apfp.new(0,0,5,-NUM_DIGITS) then
609: return Apfp.new(0,0,5,-NUM_DIGITS)
610: end
611:
612: (Apfp.new(5,-1,5,-NUM_DIGITS)* self.log_err).exp_err
613: end
# File lib/rapfp/Apfp.rb, line 52
52: def to_s
53: b = self.norm
54: if b.mant < 0 then
55: mant = -b.mant
56: sign = "-";
57: s1 = b.mant.to_s.size - 1
58: else
59: mant = b.mant
60: sign = ""
61: s1 = b.mant.to_s.size
62: end
63: if b.errmant < 0 then
64: errmant = -b.errmant
65: errsign = "-";
66: s2 = b.errmant.to_s.size - 1
67: else
68: errmant = b.errmant
69: errsign = ""
70: s2 = b.errmant.to_s.size
71: end
72:
73: sign + "0." + mant.to_s + "e" +(b.expt+s1).to_s + "+/-" + errsign + "0." + errmant.to_s + "e" + (b.errexpt+s2).to_s
74: end
# File lib/rapfp/Apfp.rb, line 337
337: def trunc
338:
339: a = self.norm
340: amant = a.mant
341: if amant == 0 then return self end
342: sz = amant.to_s.size + a.expt
343: if amant < 0 then sz_m = sz -1 else sz_m = sz end
344: if sz_m <= 0 then return @@rconst.zero end
345: sz2 = amant.size
346:
347: if sz > sz2 then (sz - sz2).times do amant *= 10 end
348: else amant = (amant.to_s[0,sz]).to_i end
349: if amant.to_s.size == 0 then amant = 0 end
350:
351: Apfp.new(amant,0,@errmant,@errexpt).norm
352: end