home *** CD-ROM | disk | FTP | other *** search
/ Chip 2002 June / Chip_2002-06_cd1.bin / zkuste / derive / download / Setup.exe / %MAINDIR% / Math / Number.mth < prev    next >
Encoding:
Text File  |  2002-05-15  |  4.9 KB  |  116 lines
  1. "File NUMBER.MTH, copyright (c) 1990-2001 by Texas Instruments Incorporated"
  2.  
  3. "The functions in this file were contributed by Johann Wiesenbauer"
  4.  
  5. "of Vienna, Austria (j.wiesenbauer@tuwien.ac.at)."
  6.  
  7. EXTENDED_GCD(a,b,q_,r_):=PROG(a:=[a,[1,0]],b:=[b,[0,1]],LOOP(IF(FIRST(b)=0,ex~
  8. it),q_:=FIRST(a)/FIRST(b),q_:=ROUND(RE(q_))+ROUND(IM(q_))*#i,r_:=a-q_*b,a:=b,~
  9. b:=r_),IF(FIRST(a)=0,RETURN([0,[0,0]])),LOOP(IF(0<=PHASE(FIRST(a))<pi/2,RETUR~
  10. N(a)),a:*#i))
  11.  
  12. SOLVE_MOD(u,x,m,a_,b_,d_):=PROG(a_:=DIF(LHS(u)-RHS(u),x),b_:=LIM(RHS(u)-LHS(u~
  13. ),x,0),d_:=GCD(a_,m),IF(MOD(b_,d_)/=0,RETURN([])),m:=m/d_,VECTOR(MOD(b_/d_*IN~
  14. VERSE_MOD(a_/d_,m),m)+k_*m,k_,0,d_-1))
  15.  
  16. CRT(a,m,m_:=1,x_:=0):=LOOP(IF(m=[],RETURN(x_)),x_:=x_+(FIRST(a)-x_)*m_*INVERS~
  17. E_MOD(m_,FIRST(m)),m_:*FIRST(m),x_:=MOD(x_,m_),a:=REST(a),m:=REST(m))
  18.  
  19. NTH_PRIME(n):=ITERATE(NEXT_PRIME(k_),k_,1,n)
  20.  
  21. PRIMEPI(x,d:=1,a:=1):=IF(d>1,SUM(IF(PRIME(d*n_+a)),n_,0,(x-a)/d),PRIMEPI(x,6,~
  22. 1)+PRIMEPI(x,6,5)+IF(x>=3,2,CHI(2,x,3,1)))
  23.  
  24. FAREY(n,a_:=1,b_:=1,l_:=[1]):=LOOP(IF(a_=b_,IF(b_=n,RETURN(SORT(l_)),[b_:+1,a~
  25. _:=1])),IF(GCD(a_,b_)=1,l_:=APPEND(l_,[a_/b_])),a_:+1)
  26.  
  27. PRIME_FACTORIZATION(n):=FACTORS(n)
  28.  
  29. DIVISORS(n):=SORT(VECTOR(PRODUCT(u_),u_,{[1]}*PRODUCT(VECTOR(MAP_LIST([v_ SUB~
  30.  1^k_],k_,{0,...,v_ SUB 2}),v_,FACTORS(n)))))
  31.  
  32. DIVISOR_SIGMA(k,n):=PRODUCT(SUM(v_ SUB 1^(k*j_),j_,0,v_ SUB 2),v_,FACTORS(n))
  33.  
  34. DIVISOR_TAU(n):=PRODUCT(v_ SUB 2+1,v_,FACTORS(n))
  35.  
  36. EULER_PHI(n):=n*PRODUCT(1-1/v_ SUB 1,v_,FACTORS(n))
  37.  
  38. PRIME_POWER?(n,k_:=1,t_:=2):=LOOP(IF(t_>n,RETURN(false)),IF(PRIME(n^(1/k_)),e~
  39. xit),k_:+1,t_:*2)
  40.  
  41. PRIMITIVE_ROOT(n,a_:=2,n_,l_,p_):=IF(n<5,n-1,IF(MOD(n,4)>0,IF(PRIME_POWER?(n/~
  42. (2-MOD(n,2))),PROG(n_:=EULER_PHI(n),l_:=(FACTORS(n_)) SUB  SUB 1,LOOP(IF(GCD(~
  43. a_,n)=1,IF(EVERY(MOD(a_^(n_/p_),n)/=1,p_,l_),RETURN(a_))),a_:+1)))))
  44.  
  45. MOEBIUS_MU(n):=PRODUCT(-MAX(2-v_ SUB 2,0),v_,FACTORS(n))
  46.  
  47. SQUAREFREE(n):=EVERY(v_ SUB 2=1,v_,FACTORS(n))
  48.  
  49. CYCLOTOMIC(n,x):=PRODUCT((x^(n/d_)-1)^MOEBIUS_MU(d_),d_,DIVISORS(n))
  50.  
  51. GEN_LUCAS(n,p,q,l0,l1):=([[0,1],[-q,p]]^n . [l0,l1]) SUB 1
  52.  
  53. U_LUCAS(n,p,q):=GEN_LUCAS(n,p,q,0,1)
  54.  
  55. V_LUCAS(n,p,q):=GEN_LUCAS(n,p,q,2,p)
  56.  
  57. U_MOD(n,p,q,m,l_,ll_,q_:=1):=PROG(m:=IF(m>0,m,0,0),l_:=MOD(2,m),ll_:=p,Output~
  58. Base:=Binary,n:=NAME_TO_CODES(STRING(n)),OutputBase:=Decimal,LOOP(IF(n=[],exi~
  59. t),IF(FIRST(n)=48,PROG(ll_:=POLY_MOD(l_*ll_-p*q_,m),l_:=POLY_MOD(l_^2-2*q_,m)~
  60. ,q_:=MOD(q_^2,m)),PROG(l_:=POLY_MOD(l_*ll_-p*q_,m),ll_:=POLY_MOD(ll_^2-2*q_*q~
  61. ,m),q_:=MOD(q*q_^2,m))),n:=REST(n)),MOD((-p*l_+2*ll_)*INVERSE_MOD(p^2-4*q,m),~
  62. m))
  63.  
  64. V_MOD(n,p,q,m,l_,ll_,q_:=1):=PROG(m:=IF(m>0,m,0,0),l_:=MOD(2,m),ll_:=p,Output~
  65. Base:=Binary,n:=NAME_TO_CODES(STRING(n)),OutputBase:=Decimal,LOOP(IF(n=[],RET~
  66. URN(l_)),IF(FIRST(n)=48,PROG(ll_:=POLY_MOD(l_*ll_-p*q_,m),l_:=POLY_MOD(l_^2-2~
  67. *q_,m),q_:=MOD(q_^2,m)),PROG(l_:=POLY_MOD(l_*ll_-p*q_,m),ll_:=POLY_MOD(ll_^2-~
  68. 2*q_*q,m),q_:=MOD(q*q_^2,m))),n:=REST(n)))
  69.  
  70. LUCAS(n):=V_LUCAS(n,1,-1)
  71.  
  72. FIBONACCI(n):=GEN_LUCAS(n,1,-1,0,1)
  73.  
  74. PELL(n):=GEN_LUCAS(n,2,-1,1,0)
  75.  
  76. LUCAS_LEHMER(p,m_):=PROG(m_:=2^p-1,SOLVE(ITERATE(MOD(s_^2-2,m_),s_,4,p-2)=0))
  77.  
  78. NEXT_MERSENNE_DEGREE(n,p_):=IF(n<2,2,PROG(p_:=NEXT_PRIME(n),LOOP(IF(LUCAS_LEH~
  79. MER(p_),RETURN(p_)),p_:=NEXT_PRIME(p_))))
  80.  
  81. MERSENNE_LIST(n):=ITERATES(NEXT_MERSENNE_DEGREE(n_),n_,2,n-1)
  82.  
  83. MERSENNE(n):=2^ITERATE(NEXT_MERSENNE_DEGREE(j_),j_,1,n)-1
  84.  
  85. MERSENNE_DEGREE(n):=[2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281~
  86. ,3217,4253,4423,9689,9941,11213,19937,21701,23209,44497,86243,110503,132049,2~
  87. 16091,756839,859433,1257787,1398269,2976221,3021377,6972593] SUB n
  88.  
  89. MERSENNE(n):=2^MERSENNE_DEGREE(n)-1
  90.  
  91. PERFECT(n):=SUBST((m_+1)*m_/2,m_,MERSENNE(n))
  92.  
  93. CONTINUED_FRACTION(u,n):=FLOOR(ITERATES(1/MOD(x_),x_,u,n))
  94.  
  95. CONVERGENT(x,k,a_,b_:=1,c_:=1,d_:=0,e_,t_):=PROG(a_:=FLOOR(x),e_:=1/MOD(x),LO~
  96. OP(IF(k=0,RETURN(a_/c_)),t_:=a_,a_:=FLOOR(e_)*a_+b_,b_:=t_,t_:=c_,c_:=FLOOR(e~
  97. _)*c_+d_,d_:=t_,e_:=1/MOD(e_),k:-1))
  98.  
  99. CONVERGENTS(x,k,a_,b_:=1,c_:=1,d_:=0,e_,l_,t_):=PROG(a_:=FLOOR(x),e_:=1/MOD(x~
  100. ),l_:=[a_],LOOP(IF(k=0,RETURN(l_)),t_:=a_,a_:=FLOOR(e_)*a_+b_,b_:=t_,t_:=c_,c~
  101. _:=FLOOR(e_)*c_+d_,d_:=t_,e_:=1/MOD(e_),k:-1,l_:=INSERT(a_/c_,l_,0)))
  102.  
  103. JACOBI(a,b,s_:=1,t_):=PROG(IF(b<0 OR MOD(b,2)=0,RETURN(?)),IF(GCD(a,b)>1,RETU~
  104. RN(0)),LOOP(IF(a=1 OR b=1,RETURN(s_)),a:=MODS(a,b),IF(a<0 AND MOD(b,4)=3,s_:=~
  105. -s_),a:=ABS(a),t_:=0,LOOP(IF(MOD(a,2)=1,exit),t_:+1,a:/2),IF(MOD(t_,2)=1 AND ~
  106. MODS(b,8)="+-"3,s_:=-s_),t_:=a,a:=b,b:=t_,IF(MOD(a,4)=3 AND MOD(b,4)=3,s_:=-s~
  107. _)))
  108.  
  109. SQUARE_ROOT(a,p,a_:=2,b_,s_:=-2,t_):=PROG(IF(p=2 OR MOD(a,p)=0,RETURN(MOD(a,p~
  110. ))),IF(JACOBI(a,p)=-1,RETURN(?)),IF(MOD(p,4)=3,RETURN(MOD(a^((p+1)/4),p))),IF~
  111. (MOD(p,8)=5,IF(MODS(a^((p-1)/4),p)=1,RETURN(MOD(a^((p+3)/8),p)),RETURN(MOD(2*~
  112. a*MOD((4*a)^((p-5)/8),p),p)))),LOOP(IF(JACOBI(a_,p)=-1,exit),a_:+1),t_:=p-1,L~
  113. OOP(s_:+1,t_:/2,IF(MOD(t_,2)=1,exit)),a_:=MOD(a_^t_,p),b_:=MOD(a^((t_+1)/2),p~
  114. ),t_:=INVERSE_MOD(a,p),LOOP(IF(s_<0,RETURN(b_)),IF(MODS((b_^2*t_)^(2^s_),p)=-~
  115. 1,b_:=MOD(a_*b_,p)),a_:=MOD(a_^2,p),s_:-1))
  116.