Sunday, 20 October 2024
Benchmarks for Vectorised Simpsons Rule
On a 16 Thread Ryzen 7. It came in slightly slower. My Global Warming Infrared Attention code, took, 2h 59m 23s non vectorised, and 3h 4mins and 42 seconds, vectorised.
Saturday, 19 October 2024
A Vectorised Simpsons Rule for JDK 23
If you are doing numerical integration in Java, try this code on the latest JDK 23, it using the
new Vectorisation methods for a faster result.
import jdk.incubator.vector.DoubleVector;
import java.util.Arrays;
import java.util.List;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicReference;
// Release Candidate Version, package will change in Release Version, true for JDK 23
import jdk.incubator.vector.DoubleVector;
import jdk.incubator.vector.VectorOperators;
import jdk.incubator.vector.VectorOperators.Operator;
import jdk.incubator.vector.VectorSpecies;
import static jdk.incubator.vector.VectorOperators.ADD;
public class SimpsonsRuleVectorised {
private static double third = 1.0/3.0;
private static int THREADS = 16;
public static double integrateThreaded(double a, double b, int N, DoubFunction f, double ...params) { // precision parameter
double h = (b - a) / (N - 1); // step size
ForkJoinPool pool = new ForkJoinPool(THREADS);
double d[] = new double[N];
for(int i=0; i {
double mul = ii%2==0? 2*third: 4*third;
if (ii==0) { mul = third; }
if (ii==N-1){ mul = third; }
double mul1 = mul;
double x = a + h * ii;
double fi = f.eval(x,ii, params);
if (Double.isNaN(fi) ){
System.err.println(f.getClass().getName() + "IS NaN at "+x);
}
d[ii] = mul1*fi;
});
}
try {
pool.shutdown();
if (!pool.awaitTermination(1, TimeUnit.HOURS)){
pool.shutdownNow();
}
} catch (Exception e){}
DoubleVector doubleVector = DoubleVector.fromArray(DoubleVector.SPECIES_64, d, 0);
double sum = doubleVector.reduceLanes(ADD );
return sum * h;
}
public static double integrate(double a, double b, int N, DoubFunction f, double ...params) { // precision parameter
double h = (b - a) / (N - 1); // step size
double fa = f.eval(a,0, params);
double fb = f.eval(b, N-1, params);
if (Double.isNaN(fa) ){
System.err.println(f.getClass().getName() + "IS NaN at "+a);
}
if (Double.isNaN(fb)){
System.err.println(f.getClass().getName() + "IS NaN at "+b);
}
// 1/3 terms
double sum = third * (fa + fb);
// 4/3 terms
for (int i = 1; i < N - 1; i += 2) {
double x = a + h * i;
double fx = f.eval(x,i, params);
if (Double.isNaN(fx)){
System.err.println(f.getClass().getName() + "IS NaN at "+x);
}
sum += 4.0 * third * fx;
}
// 2/3 terms
for (int i = 2; i < N - 1; i += 2) {
double x = a + h * i;
double fx = f.eval(x,i, params);
if (Double.isNaN(fx)){
System.err.println(f.getClass().getName() + "IS NaN at "+x);
}
sum += 2.0 * third * fx;
}
return sum * h;
}
public static double integrateConsecutive(double a, double b, int N, DoubFunction f, double ...params) { // precision parameter
double h = (b - a) / (N - 1); // step size
double fa = f.eval(a,0, params);
if (Double.isNaN(fa) ){
System.err.println(f.getClass().getName() + "IS NaN at "+a);
}
// 1/3 terms
double sum = third*fa;
// 4/3 terms
boolean isOdd = true;
for (int i = 1; i < N - 1; i += 1) {
double x = a + h * i;
double fx = f.eval(x,i, params);
if (Double.isNaN(fx)){
System.err.println(f.getClass().getName() + "IS NaN at "+x);
}
if (isOdd) {
sum += 4.0 * third * fx;
} else {
sum += 2.0 * third * fx;
}
isOdd=!isOdd;
}
double fb = f.eval(b, N-1, params);
if (Double.isNaN(fb)){
System.err.println(f.getClass().getName() + "IS NaN at "+b);
}
sum = sum+ fb*third;
return sum * h;
}
public static void main(String args[]){
// Roots of polynumerial to integrate
List in = Arrays.asList(-0.9, -0.8,-0.7, -0.6,-0.5, -0.4, -0.3, -0.2, -0.1, 0, .1,.2,.3,.4,.5, .6,.7, .8, .9 );
DoubFunction func = new DoubFunction() {
@Override
double evalInner(double x, double[] params, int i) {
return in.stream().map(y->y.doubleValue()-x).reduce(1.0,(a,b)->(a*b));
}
};
double consec=0;
long startConsec = System.currentTimeMillis();
for(int i=1;i<1000; i++) {
consec = integrateConsecutive(-1, 1, 100000, func);
}
double timeConsec = (System.currentTimeMillis() - startConsec)/1000.0;
double standard=0;
long startStandard = System.currentTimeMillis();
for(int i=1; i<1000; i++) {
standard = integrate(-1,1,100000,func);
}
double timeStandard = (System.currentTimeMillis() - startStandard)/1000.0;
double threaded=0;
long startThreaded = System.currentTimeMillis();
for(int i=1;i<1000; i++) {
threaded = integrateThreaded(-1, 1, 100000, func);
}
double timeThreaded = (System.currentTimeMillis() - startThreaded)/1000.0;
System.out.println("Standard Integrator: "+standard+" time taken: "+timeStandard+" seconds");
System.out.println("Consecutive Integrator: "+consec+" time taken: "+timeConsec+" seconds");
System.out.println("threaded Integrator: "+threaded+" time taken: "+timeThreaded+" seconds");
}
}
class AtomicDoubleLocal77 {
private AtomicReference value = new AtomicReference(Double.valueOf(0.0));
double getAndAdd(double delta) {
while (true) {
Double currentValue = value.get();
Double newValue = Double.valueOf(currentValue.doubleValue() + delta);
if (value.compareAndSet(currentValue, newValue))
return currentValue.doubleValue();
}
}
}
Thursday, 17 October 2024
Recent Paper has new limit on neutrino - majoron interactions from the SN1987 Supernova
The supernova of 1987, SN1987a which occurred in the large magilangic cloud, is the only supernova where neutrinos (some 24), interactions on earth have been recorded. These recorded interactions have
now allow the researcher, P.I. Ballesteros and Christa Volpe, https://arxiv.org/abs/2410.11517, to limit potential interactions between neutrinos and majorons (a Majoron is a spin zero particle like a light Huggs, giving mass just to the neutrino).
The limits are force contant of around 10^-7 compare. This might also limit an axial force, although that would be a spin-1 psuedovector particle, previously we look at a force contant in the range a few*10^-5 so this
paper might simiilar reduce the limits on strength of how force by a factor of a hundred. The previous work on SN1987, from the year 2000, and published in Phy Rev D, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.62.023004
limited the majoron interactiom in the range, 3×10−7≲𝑔≲2×10−5 or 𝑔≳3×10−4 so left open the range 2*10^-5 to 3*10-4 allow or orignal guess of force constant (making the weak assumption that majoron and axi-photon limits would be similar.)
Monday, 23 September 2024
NSI to solve the KOTO anomaly
This article in 2008, shows an attempted to resolve the KOTO and Invisible beutry decay anomalies via a new neutrino interaction, (there U(1) B-L)
https://arxiv.org/abs/2008.09793
Differences between Nova and T2K experiment hint (1.8 sigma) at a Neutrino Firth Force
https://arxiv.org/pdf/2409.10599 states the current status on a NSI is ~1.8 sigma and it is due to a CP violating phase Tension in measurements.
With the axial force C is violated in matter neutrino interactions, due to matter being matter and not anti-matter, P and perhaps CP might
be violated if the matter has a net axial charge where the neutrinos where passing.
Wednesday, 28 August 2024
MOre invisible beauty decays.
https://arxiv.org/abs/2312.12507 In a recent update to the invisable beauty decay found at Bella, finds it favours a two body decay, (but with final mass 0.6 GeV), our
our axial force has 6 invisable end products, electron, muon and tau neutrino anti-neutrino pairs, but also there near sterile right handed complements, could the right handed tau neutrino have
a 600MeV mass? If so it would make a strong decay matter candidate, if and only if its does't decay or its decay is very slow.
Tuesday, 28 November 2023
Subscribe to:
Posts (Atom)