/usr/share/yacas/numbers.rep/NumberTheory.ys is in yacas 1.3.3-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 | /* Implementation of some number theoretical functions for Yacas */
/* (C) 2002 Pablo De Napoli <pdenapo@yahoo.com> under GNU GPL */
/* DivisorsList(n) = the list of divisors of n */
DivisorsList(n_IsPositiveInteger) <--
[
Local(nFactors,f,result,oldresult,x);
nFactors:= Factors(n);
result := {1};
ForEach (f,nFactors)
[
oldresult := result;
For (k:=1,k<=f[2],k++)
ForEach (x,oldresult)
result:=Append(result,x*f[1]^k);
];
result;
];
/* This function performs a sum where sumvar runs through
the divisors of n
For example SumForDivisors(d,10,d^2)
sums d^2 with d walking through the divisors of 10
LocalSymbols is needed since we use Eval() inside
Look at Programming in Yacas: Evaluating Variables in the Wrong
Scope */
Function ("SumForDivisors",{sumvar,n,sumbody}) LocalSymbols(s,d)
[
Local(s,d);
s:=0;
ForEach (d,DivisorsList(n))
[
MacroLocal(sumvar);
MacroSet(sumvar,d);
s:=s+Eval(sumbody);
];
s;
];
UnFence("SumForDivisors",3);
HoldArg("SumForDivisors",sumvar);
HoldArg("SumForDivisors",sumbody);
/* Returns a list of the square-free divisors of n */
SquareFreeDivisorsList(n_IsPositiveInteger) <--
[
Local(nFactors,f,result,oldresult,x);
nFactors:= Factors(n);
result := {1};
ForEach (f,nFactors)
[
oldresult := result;
ForEach (x,oldresult)
result:=Append(result,x*f[1]);
];
result;
];
/* Returns a list of pairs {d,m}
where d runs through the square free divisors of n
and m=Moebius(m)
This is much more efficient than making a list of all
square-free divisors of n, and then compute Moebius on each of them.
It is useful for computing the Cyclotomic polinomials.
It can be useful in other computations based on
Moebius inversion formula. */
MoebiusDivisorsList(n_IsPositiveInteger) <--
[
Local(nFactors,f,result,oldresult,x);
nFactors:= Factors(n);
result := {{1,1}};
ForEach (f,nFactors)
[
oldresult := result;
ForEach (x,oldresult)
result:=Append(result,{x[1]*f[1],-x[2]});
];
result;
];
/* RamanujanSum(k,n) = the sum of the n-th powers of the
k-th primitive roots of the identity */
10 # RamanujanSum(k_IsPositiveInteger,0) <-- Totient(k);
20 # RamanujanSum(k_IsPositiveInteger,n_IsPositiveInteger) <--
[
Local(s,gcd,d);
s:= 0;
gcd := Gcd(n,k);
ForEach (d,DivisorsList(gcd))
s:=s+d*Moebius(k/d);
s;
];
/** Compute the Jacobi symbol JS(m/n) - n must be odd, both positive.
See the Algo book for documentation.
*/
10 # JacobiSymbol(_a, 1) <-- 1;
15 # JacobiSymbol(0, _b) <-- 0;
18 # JacobiSymbol(_a, _b) _ (Gcd(a,b)>1) <-- 0;
20 # JacobiSymbol(_a, b_IsOdd)_(a>=Abs(b) Or a<0) <-- JacobiSymbol(Mod(a,Abs(b)),Abs(b));
30 # JacobiSymbol(a_IsEven, b_IsOdd) <--
[
Local(c, s);
// compute c,s where a=c*2^s and c is odd
{c,s}:=FindPrimeFactorSimple(a, 2); // use the "Simple" function because we don't expect a worst case here
If(Mod(s,2)=1 And Abs(Mod(b,8)-4)=1, -1, 1) * JacobiSymbol(c,b);
];
40 # JacobiSymbol(a_IsOdd, b_IsOdd) <-- If(Mod(a,4)=3 And Mod(b,4)=3, -1, 1) * JacobiSymbol(b,a);
|