Factorial

Swatch Watch · Post by **Swatch Watch** » Sun Jul 19, 2009 9:06 am

The lower code works ok but if N is e.g 150 it gives back a wrong factorial. So my question is: Do you know any more correct and or faster way to get a factorial? I need these high factorials for the Bernstein polynomial.

Code: Select all

Procedure.d Factorial(N.b)
E.d=1
For S=2 To N
E*S
Next
ProcedureReturn E
EndProcedure

N = 120
MessageRequester("Factorial of "+Str(N),StrD(Factorial(N),0))

Bernstein polynomial:

Code: Select all

Procedure.d Factorial(N.l)
E.d=1
For S=2 To N
E*S
Next
ProcedureReturn E
EndProcedure

Global NewList BernsteinElemente()

Procedure Bernstein(Width,Height,FrontColor,BackColor)
ResetList(BernsteinElemente())
NextElement(BernsteinElemente())
N = BernsteinElemente()

ImageID = CreateImage(#PB_Any,Width,Height)
StartDrawing(ImageOutput(ImageID))
Box(0,0,Width,Height,BackColor)
For NV = 1 To ListSize(BernsteinElemente())-1
NextElement(BernsteinElemente())
V = BernsteinElemente()
Old = Height
For X = 1 To Width
New = Height-Height*(Factorial(N)/(Factorial(V)*Factorial(N-V)))*Pow((X/Width),V)*Pow(1-(X/Width),N-V)
LineXY(X-1,Old,X,New,FrontColor)
Old = New
Next
Next
StopDrawing()
ProcedureReturn ImageID
EndProcedure


N = 30
AddElement(BernsteinElemente()) : BernsteinElemente() = N
For W  = 0 To N
AddElement(BernsteinElemente()) : BernsteinElemente() = W
Next
SetClipboardImage(Bernstein(500,500,RGB(0,120,220),RGB(20,20,20)))

Kaeru Gaman · Post by **Kaeru Gaman** » Sun Jul 19, 2009 9:41 am

you can't use the standard floating point storage for such high numbers.

the mantissa of a double has 52bit, so the edge of precision is 2^53, not more.
http://en.wikipedia.org/wiki/Floating_p ... esentation
if you store higher numbers, the rest is filled with zeros when displayed.

you would need some high-number-library to make your calculations, or take a totally different approach.

dioxin · Post by **dioxin** » Sun Jul 19, 2009 5:38 pm

The FPU can handle 80bit Extended floating point numbers which will give you up to 1754! before it overflows but you'll have to program it in ASM as PB doesn't handle extended FP values. Psuedo code follows.

Code: Select all

!fld1           
!fld1
!fld Argument        'the input argument

lp:
!fmul st(1),st(0)
!fsub st(0),st(2)
!fcomi st(0),st(2)
!jne lp
!fstp st(0)
!fstp answer         'where to store the result, argument!
!fstp st(0)

You would probably choose not to store the result (as PB has no variable type to handle it) but to continue the calculation in the FPU using ASM.

You can also do the calculation using logs and this will allow you to calculate extremely large factorials:

Code: Select all

FOR r = 1 TO argument
    y = y + LOG10(r)
NEXT
'y now = log10(argument!)

Depending on the values of N and V in your code you may also be able to avoid the overflow by noting that
n!/(v! * (n-v)!) = n!/v! / (n-v)! and n!/v! can be calculated by getting the product of integers from v+1 to n instead of 1 to n.

So instead of:

Code: Select all

For S=2 To N 
E*S 
Next

You'd do:

Code: Select all

For S=V+1 To N 
E*S 
Next

which will result in a smaller product.

Kaeru Gaman · Post by **Kaeru Gaman** » Sun Jul 19, 2009 9:11 pm

@dioxin

the Problem is not the overflow, but the precision.

an 80bit Float may hold 1754!, but not to the presision on the last digit.

120! sure can be kept in a Double, a 64bit Float, but the precision is lost.

Swatch Watch · Post by **Swatch Watch** » Mon Jul 20, 2009 9:39 pm

Thank you for the answers.

Frarth · Post by **Frarth** » Tue Jul 21, 2009 11:23 am

I have the same problem converting a Double from string to value:

result.d = ValD("0.00057")

returns 0.00056999999...

Why is the precision lost and is there a way to avoid this?

Frank

Helle · Post by **Helle** » Tue Jul 21, 2009 11:46 am

This is a little precision-check for 80-bit-factorials:

Code: Select all

;Precision-Test Factorial 80-Bit

Global Mantissa.d
Global Test.d
Global Exponent.l
Global Argument.l = 1754     ;Maximum for 80-Bit, Range 2-1754

!fninit                 ;64-bit precision 
;calculate the factorial
!fld1            
!fld1
!fld1
!fadd st0,st1
!mov ecx,[v_Argument]       
!dec ecx
!@@: 
!fmul st1,st0
!fadd st0,st2
!dec ecx
!jnz @b
!fstp st0  
;conversion in decimal-exponents-output, based of code from akj
!fxtract
!fxch
!fld st0
!fldlg2
!fmulp
!frndint
!fist [v_Exponent]  
!fldl2t  
!fmulp 
!fsubp  
!fst qword[v_Test]
If Abs(Test) < 1
  !f2xm1
  !fld1
  !faddp
 Else
  !fld1 
  !fld st0
  !faddp
  !fdivp
  !f2xm1
  !fld1
  !faddp
  !fld st0
  !fmulp
EndIf 
!fmulp   
!fstp qword[v_Mantissa]

!fstp st0               ;all st free

While Mantissa < 1.0 : Mantissa * 10.0 : Exponent - 1 : Wend 
While Mantissa >= 10.0 : Mantissa / 10.0 : Exponent + 1 : Wend 

Factorial$ = "80-Bit-Test PB            " + StrD(Mantissa, 20) + "e+" + Str(Exponent) 
WinCalc$ = "Windows-Calculator : 1.9792618901050100553817943275326e+4930" 

MessageRequester("Precision-Test Factorial 80-Bit for Value 1754", WinCalc$ + #LFCR$ + Factorial$)

I think, the precision is o.K. for the most purposes!

Gruss
Helle

gnasen · Post by **gnasen** » Tue Jul 21, 2009 11:56 am

Frarth wrote:I have the same problem converting a Double from string to value:

result.d = ValD("0.00057")

returns 0.00056999999...

Why is the precision lost and is there a way to avoid this?

Frank

The problem is, that a float (double etc) number stores the value in a special way.

The number itself is turned into a binary number, 24 goes to 11000.
If you now try to store 2,4 it makes the following:
It stores 11000 and adds the information: Shift the comma one diggit to the right: 24 -> 2,4

This is very good to store even very big or small numbers.
However there is some problem with the precision.
First, the number and the shift are limited by the variable size.
Second, you cant display every number with finite diggits of a binary number.
You can try it yourself to convert 10/3 to a binary.
The result wil be the infinite number
11.010101010101...
so weve got some problem here.

I hope it helped you a little bit

Kaeru Gaman · Post by **Kaeru Gaman** » Tue Jul 21, 2009 12:36 pm

is there a way to avoid this?

only by using fixcomma, this means storing an Integer and calculating with e.g. 1/1000

you can store a 12.00057 as a 1200057, and add the comma while displaying.

disadvantage: fixcomma numbers are strictly limited to the precision you define them.
when you use factor 100000 as above, the sixth digit never can be stored.
you can store numbers from -21 474.83648 to +21 474.83647 in a Long
and numbers from -92 233 720 368 547.75808 to +92 233 720 368 547.75807 in a Quad

Arctic Fox · Post by **Arctic Fox** » Tue Jul 21, 2009 4:59 pm

APM-libraries like pbAPM or MAPM can handle very big numbers without loss of precision.

pbAPM: http://purearea.net/pb/download/userlibs/pbAPM.zip
MAPM: http://www.purebasic.fr/english/viewtopic.php?t=35173 (thanks again jack!

)

Note that pbAPM is not working with the latest versions of PureBasic, though

Frarth · Post by **Frarth** » Tue Jul 21, 2009 6:56 pm

Thankx everyone! Yes, this helps a lot.

Frank