Your Ad Here Your Ad Here
Results 1 to 7 of 7

Thread: Number Scrunching

  1. #1
    MediaSlayer's Avatar slowly going deaf
    Join Date
    May 2003
    Location
    ur anus
    Posts
    2,783
    what does it mean


    sending fiery missiles in manker's japan's general direction.

  2. Lounge   -   #2
    LeGoMyFnLeg
    Guest
    Number crunching.

    One of the most famous examples.
      program piqp

    c  This program employs the recently discovered digit extraction scheme
    c  to produce hex digits of pi.  This code is valid up to ic = 2^24 on
    c  systems with IEEE arithmetic.

          real*8 pid, s1, s2, s3, s4, series
          parameter (ic = 100000, nhx = 12)
          character*1 chx(nhx)

    c  ic is the hex digit position -- output begins at position ic + 1.

          s1 = series (1, ic)
          s2 = series (4, ic)
          s3 = series (5, ic)
          s4 = series (6, ic)
          pid = mod (4.d0 * s1 - 2.d0 * s2 - s3 - s4, 1.d0) + 1.d0
          call ihex (pid, nhx, chx)
          write (6, 1) ic, pid, chx
    1    format ('Pi hex digit computation'/'position =',i12,' + 1'/
        $  f20.15/12a1)

          stop
          end

          subroutine ihex (x, nhx, chx)

    c  This returns, in chx, the first nhx hex digits of the fraction of x.

          real*8 x, y
          character*1 chx(nhx), hx(0:15)
          data hx/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A',
        $  'B', 'C', 'D', 'E', 'F'/

          y = abs (x)

          do i = 1, nhx
            y = 16.d0 * (y - aint (y))
            chx(i) = hx(int(y))
          enddo

          return
          end

          function series (m, ic)

    c  This routine evaluates the series  sum_k 16^(ic-k)/(8*k+m)
    c  using the modular exponentiation technique.

          real*8 ak, eps, expm, p, s, series, t
          parameter (eps = 1d-17)

          s = 0.d0

    c  Sum the series up to ic.

          do k = 0, ic - 1
            ak = 8 * k + m
            p = ic - k
            t = expm (p, ak)
            s = mod (s + t / ak, 1.d0)
          enddo

    c  Compute a few terms where k >= ic.

          do k = ic, ic + 100
            ak = 8 * k + m
            t = 16.d0 ** (ic - k) / ak
            if (t .lt. eps) goto 100
            s = mod (s + t, 1.d0)
          enddo

    100  continue

          series = s
          return
          end

          function expm (p, ak)

    c  expm = 16^p mod ak.  This routine uses the left-to-right binary
    c&nbsp; exponentiation scheme.&nbsp; It is valid for&nbsp; ak <= 2^24.

    &nbsp; &nbsp; &nbsp; parameter (ntp = 25)
    &nbsp; &nbsp; &nbsp; real*8 ak, expm, p, p1, pt, r, tp(ntp)
    &nbsp; &nbsp; &nbsp; save tp
    &nbsp; &nbsp; &nbsp; data tp/ntp*0.d0/

    c&nbsp; If this is the first call to expm, fill the power of two table tp.

    &nbsp; &nbsp; &nbsp; if (tp(1) .eq. 0.d0) then
    &nbsp; &nbsp; &nbsp; &nbsp; tp(1) = 1.d0

    &nbsp; &nbsp; &nbsp; &nbsp; do i = 2, ntp
    &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; tp(i) = 2.d0 * tp(i-1)
    &nbsp; &nbsp; &nbsp; &nbsp; enddo
    &nbsp; &nbsp; &nbsp; endif

    &nbsp; &nbsp; &nbsp; if (ak .eq. 1.d0) then
    &nbsp; &nbsp; &nbsp; &nbsp; expm = 0.d0
    &nbsp; &nbsp; &nbsp; &nbsp; return
    &nbsp; &nbsp; &nbsp; endif

    c&nbsp; Find the greatest power of two less than or equal to p.

    &nbsp; &nbsp; &nbsp; do i = 1, ntp
    &nbsp; &nbsp; &nbsp; &nbsp; if (tp(i) .gt. p) goto 100
    &nbsp; &nbsp; &nbsp; enddo

    100&nbsp; pt = tp(i-1)
    &nbsp; &nbsp; &nbsp; p1 = p
    &nbsp; &nbsp; &nbsp; r = 1.d0

    c&nbsp; Perform binary exponentiation algorithm modulo ak.

    110&nbsp; continue

    &nbsp; &nbsp; &nbsp; if (p1 .ge. pt) then
    &nbsp; &nbsp; &nbsp; &nbsp; r = mod (16.d0 * r, ak)
    &nbsp; &nbsp; &nbsp; &nbsp; p1 = p1 - pt
    &nbsp; &nbsp; &nbsp; endif
    &nbsp; &nbsp; &nbsp; pt = 0.5d0 * pt
    &nbsp; &nbsp; &nbsp; if (pt .ge. 1.d0) then
    &nbsp; &nbsp; &nbsp; &nbsp; r = mod (r * r, ak)
    &nbsp; &nbsp; &nbsp; &nbsp; goto 110
    &nbsp; &nbsp; &nbsp; endif

    &nbsp; &nbsp; &nbsp; expm = mod (r, ak)

    &nbsp; &nbsp; &nbsp; return
    &nbsp; &nbsp; &nbsp; end

    Another example

    http://www.activemind.com/Mysterious/Topic...e_equation.html

  3. Lounge   -   #3
    Poster
    Join Date
    Aug 2003
    Posts
    1,477
    You could have just used a dictionary defenition...

  4. Lounge   -   #4
    LeGoMyFnLeg
    Guest
    Feel free to add to the discussion, I believe that&#39;s what forums are for?

  5. Lounge   -   #5
    Poster
    Join Date
    Aug 2003
    Posts
    1,477
    Nah I prefer to drop a little comment here and there and then walk off again.

  6. Lounge   -   #6
    LeGoMyFnLeg
    Guest
    Should I provide a definition for that?

  7. Lounge   -   #7
    Poster
    Join Date
    Aug 2003
    Posts
    1,477

Bookmarks

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •