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

Methods

*   **   +   -   /   <   <=   ==   >   >=   about_log10   abs   acos   asin   asin_basic   atan   atan2   cos   cosh   display_val   eq   erf   errexpt   errexptadd!   errmant   errmanttimes10!   exp   exp_err   expt   exptadd!   frac   ge   gt   le   log   log10   log10_err   log_basic   log_basic_err   log_basic_init   log_err   lt   makerr   mant   manttimes10!   neg   new   norm   period   seterr!   setrconst   sin   sinh   sqrt   sqrt_err   tan   tanh   to_f   to_s   trunc  

Public Class methods

@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.

[Source]

    # 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

Public Instance methods

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # File lib/rapfp/Apfp.rb, line 229
229: def -(other)
230:   (self+other.neg)
231: end

[Source]

     # 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

[Source]

     # File lib/rapfp/Apfp.rb, line 292
292: def <(other)
293: 
294:   (self - other).mant < 0
295: 
296: end

[Source]

     # File lib/rapfp/Apfp.rb, line 314
314: def <=(other)
315: 
316:   !(self > (other))
317: 
318: end

[Source]

     # File lib/rapfp/Apfp.rb, line 304
304: def ==(other)
305:   (self - other).mant == 0
306: end

[Source]

     # File lib/rapfp/Apfp.rb, line 298
298: def >(other)
299: 
300:   other < self
301: 
302: end

[Source]

     # File lib/rapfp/Apfp.rb, line 308
308: def >=(other)
309: 
310:   !(self < (other))
311: 
312: end

very rough & fast

[Source]

     # File lib/rapfp/Apfp.rb, line 524
524: def about_log10
525:   b = self.norm
526:   return ap_int(b.expt)
527: end

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # File lib/rapfp/Apfp.rb, line 326
326: def eq(other)
327:   !(self < other) && !(self > other)
328: end

[Source]

     # 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

[Source]

    # File lib/rapfp/Apfp.rb, line 89
89: def errexpt
90:   @errexpt
91: end

[Source]

     # File lib/rapfp/Apfp.rb, line 165
165: def errexptadd!(other)
166:   @errexpt += other
167:   self
168: end

[Source]

    # File lib/rapfp/Apfp.rb, line 86
86: def errmant
87:   @errmant
88: end

[Source]

     # File lib/rapfp/Apfp.rb, line 155
155: def errmanttimes10!
156:   @errmant *= 10
157:   self
158: end

[Source]

     # 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

[Source]

     # 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

[Source]

    # File lib/rapfp/Apfp.rb, line 83
83: def expt
84:   @expt
85: end

[Source]

     # File lib/rapfp/Apfp.rb, line 160
160: def exptadd!(other)
161:   @expt += other
162:   self
163: end

[Source]

     # File lib/rapfp/Apfp.rb, line 354
354: def frac
355: 
356:   (self - self.trunc).norm
357: 
358: end

[Source]

     # File lib/rapfp/Apfp.rb, line 332
332: def ge(other)
333:   !(self < other)
334: end

[Source]

     # File lib/rapfp/Apfp.rb, line 323
323: def gt(other)
324:   (self-self.makerr) < (other+other.makerr)
325: end

[Source]

     # File lib/rapfp/Apfp.rb, line 329
329: def le(other)
330:   !(self > other)
331: end

log base e of any positive value

[Source]

     # File lib/rapfp/Apfp.rb, line 559
559: def log
560: #  puts "log e of " + self.to_s
561: 
562:   b = self.log10*@@rconst.log_e_10
563: #  puts "log(" + self.to_s + "} = " + b.to_s
564:   return b
565: 
566: end

log base 10 of positive value

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

log base e of positive value for purpose of error

[Source]

     # File lib/rapfp/Apfp.rb, line 569
569: def log_err
570:   return self.log10_err*Apfp.new(2302585093,-9,5,-9)
571: end

[Source]

     # File lib/rapfp/Apfp.rb, line 320
320: def lt(other)
321:   (self+self.makerr) < (other-other.makerr)
322: end

[Source]

     # File lib/rapfp/Apfp.rb, line 256
256: def makerr
257:   Apfp.new(@errmant,@errexpt,5,-NUM_DIGITS).norm
258: end

[Source]

    # File lib/rapfp/Apfp.rb, line 80
80: def mant
81:   @mant
82: end

[Source]

     # File lib/rapfp/Apfp.rb, line 151
151: def manttimes10!
152:   @mant *= 10
153:   self
154: end

[Source]

     # File lib/rapfp/Apfp.rb, line 193
193: def neg
194:   mant = -1*@mant
195:   Apfp.new(mant,@expt,@errmant.abs,@errexpt)
196: end

[Source]

     # 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

[Source]

     # 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

[Source]

     # File lib/rapfp/Apfp.rb, line 260
260: def seterr!(limit)
261:   @errmant = limit.mant
262:   @errexpt = limit.expt
263:   self.norm
264: end

[Source]

    # File lib/rapfp/Apfp.rb, line 32
32: def setrconst(c)
33:   @@rconst = c
34: end

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # 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

[Source]

     # File lib/rapfp/Apfp.rb, line 575
575: def tan
576: 
577:   self.sin/self.cos
578: 
579: end

[Source]

     # File lib/rapfp/Apfp.rb, line 638
638: def tanh
639: 
640:         self.sinh/self.cosh
641: 
642: end

[Source]

    # File lib/rapfp/Apfp.rb, line 76
76: def to_f
77:   self.to_s.to_f
78: end

[Source]

    # 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

[Source]

     # 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

[Validate]