flat assembler
Message board for the users of flat assembler.

Index > Projects and Ideas > Hex digit extraction of PI:

Author
Thread Post new topic Reply to topic
bitRAKE



Joined: 21 Jul 2003
Posts: 4046
Location: vpcmpistri
bitRAKE 11 Nov 2007, 06:32
I just coded and verified (10 000 digits, with two other sources) this algorithm. It isn't good for anything, but maybe someone else might like to play with it.

It is a direct conversion of the BBP algorithm. Here are some single digit timings on my Pentium M:
Code:
      digit   time
----------------------------
        10^4      16 ms
     10^5     219 ms
     10^6    2860 ms
     10^7   33938 ms
     10^8  407828 ms    
An SSE2 version might be interesting to code in the future...
Code:
; 16^EBP MOD ECX=8k+z
; Special cases:
;       EBP=0, return 1
;    ECX=1, return 0
;
; EDX is return value
mpow:      xor edx,edx
 cmp ecx,1
   je .x
       bsr ebx,ebp
 mov edx,1
   jne .1
.x:       retn

.0:     mul eax
     div ecx
.1:      bt ebp,ebx
  mov eax,16
  jnc .2
      mul edx
     div ecx
.2:      dec ebx
     mov eax, edx
        jns .0
      retn

series: xor edi,edi     ; fraction
.1:   call mpow       ; 16^EBP mod ECX
    xor eax,eax
 div ecx         ; EDX:0 / 8k+m
  add ecx,8       ; 8k+m, k[0,EBP]
    add edi,eax
 dec ebp
     jns .1

  mov ebp,10000000h
.2:    mov eax,ebp
 xor edx,edx
 div ecx
     add ecx,8
   shr ebp,4
   cmp ecx,ebp
 lea edi,[edi+eax]
   jna .2
      retn

PiHexDigit:
 pushad
      mov ebp,[esp+36]
    mov ecx,1
   call series
 lea esi,[edi*4]

 mov ebp,[esp+36]
    mov ecx,4
   call series
 sub esi,edi
 sub esi,edi

     mov ebp,[esp+36]
    mov ecx,5
   call series
 sub esi,edi

     mov ebp,[esp+36]
    mov ecx,6
   call series
 sub esi,edi
 mov [esp+28],esi
    popad
       retn 4

; use it like this...
 invoke PiHexDigit, 1000000
; top nibble of EAX is hex digit    
[3.]243F6 A8885 A308D 31319 8A2E0 37073 44A40 93822 299F3 1D008 2EFA9 8EC4E 6C894 52821 E638D 01377 ...

Typically, each call returns at least five valid hex digits, but I haven't tested the validity of that claim rigorously.
Post 11 Nov 2007, 06:32
View user's profile Send private message Visit poster's website Reply with quote
edfed



Joined: 20 Feb 2006
Posts: 4347
Location: Now
edfed 13 Nov 2007, 14:20
good job!
Post 13 Nov 2007, 14:20
View user's profile Send private message Visit poster's website Reply with quote
bitRAKE



Joined: 21 Jul 2003
Posts: 4046
Location: vpcmpistri
bitRAKE 13 Nov 2007, 16:52
Next step is to remove the divisions with Montgomery exponentiation. Then I'll re-examine what is required for an SSE2 implementation. Rounding is free with scaled integers, so I doubt moving to floating point would help.
Post 13 Nov 2007, 16:52
View user's profile Send private message Visit poster's website Reply with quote
Display posts from previous:
Post new topic Reply to topic

Jump to:  


< Last Thread | Next Thread >
Forum Rules:
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
You cannot attach files in this forum
You can download files in this forum


Copyright © 1999-2024, Tomasz Grysztar. Also on GitHub, YouTube.

Website powered by rwasa.